diff --git a/.circleci/config.yml b/.circleci/config.yml index 9ebfa314..00154863 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -37,7 +37,7 @@ jobs: name: Run tests # This assumes pytest is installed via the install-package step above command: pytest - + # Invoke jobs via workflows # See: https://circleci.com/docs/2.0/configuration-reference/#workflows workflows: diff --git a/.github/workflows/automated-testing.yml b/.github/workflows/automated-testing.yml index 510a73f5..48eac0dd 100644 --- a/.github/workflows/automated-testing.yml +++ b/.github/workflows/automated-testing.yml @@ -1,33 +1,41 @@ -# This workflow will install Python dependencies and run tests with a single version of Python -# For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions - -name: Python application - on: push: - branches: [ "main" ] + branches: + - main pull_request: - branches: [ "main" ] - -permissions: - contents: read - + branches: + - main +name: CI jobs: build: - runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - name: Set up Python 3.10 + uses: actions/setup-python@v3 + with: + python-version: "3.10" + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install pytest + python -m pip install . + - name: Test with pytest + run: | + python -m pytest + flake8-lint: + name: flake8 + runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 - - name: Set up Python 3.10 - uses: actions/setup-python@v3 - with: - python-version: "3.10" - - name: Install dependencies - run: | - python -m pip install --upgrade pip - pip install flake8 pytest - python -m pip install . - - name: Test with pytest - run: | - python -m pytest + - uses: actions/checkout@v3 + - name: Set up Python environment + uses: actions/setup-python@v4 + with: + python-version: "3.10" + # Install flake8 extensions (this step is not required. Default is "None"). + - name: Set up flake8 + run: pip install flake8-quotes + - name: flake8 Lint + run: | + flake8 hogben diff --git a/hogben/__init__.py b/hogben/__init__.py index fcf90872..cbff3dcb 100644 --- a/hogben/__init__.py +++ b/hogben/__init__.py @@ -1,3 +1,3 @@ """Top-level package for HOGBEN""" -__version__ = "1.2.1" +__version__ = '1.2.1' diff --git a/hogben/angles.py b/hogben/angles.py index 701daff7..fe2758b7 100644 --- a/hogben/angles.py +++ b/hogben/angles.py @@ -1,5 +1,4 @@ import os -import sys import time import numpy as np @@ -16,18 +15,14 @@ def _angle_results_visualise(save_path): save_path (str): path to directory to save results to. """ - from models.samples import similar_sld_sample_1, similar_sld_sample_2 - from models.samples import thin_layer_sample_1, thin_layer_sample_2 - from models.samples import simple_sample, many_param_sample - from models.bilayers import BilayerDMPC, BilayerDPPC - from models.magnetic import SampleYIG + from models.samples import simple_sample # Choose sample here. sample = simple_sample() # Choose contrasts here (only bilayers should use this). contrasts = [] - #contrasts = [-0.56, 6.36] + # contrasts = [-0.56, 6.36] # Number of points and counting time for the initial angle choice. points = 100 @@ -37,13 +32,13 @@ def _angle_results_visualise(save_path): angle_range = np.linspace(0.2, 4, 500) initial_angle = angle_choice(sample, [], angle_range, points, time, save_path, 'initial', contrasts) - # Plot how the choice of angle changes with initial angle counting time. angle_range = np.linspace(0.2, 4, 50) - time_range = np.linspace(0, time*8, 50) + time_range = np.linspace(0, time * 8, 50) angle_choice_with_time(sample, initial_angle, angle_range, time_range, points, time, save_path, contrasts) + def _angle_results_optimise(save_path): """Optimises the choice measurement angles and counting times for a sample. @@ -51,21 +46,16 @@ def _angle_results_optimise(save_path): save_path (str): path to directory to save results to. """ - from models.samples import similar_sld_sample_1, similar_sld_sample_2 - from models.samples import thin_layer_sample_1, thin_layer_sample_2 - from models.samples import simple_sample, many_param_sample - from models.bilayers import BilayerDMPC, BilayerDPPC - from models.magnetic import SampleYIG - + from models.samples import simple_sample # Choose sample here. sample = simple_sample() # Choose contrasts here (only bilayers should use this). contrasts = [] - #contrasts = [-0.56, 6.36] + # contrasts = [-0.56, 6.36] # Total time budget. - total_time = 1000 # A large time improves DE convergence. + total_time = 1000 # A large time improves DE convergence. # Interval containing angles to consider. angle_bounds = (0.2, 4.0) @@ -73,7 +63,7 @@ def _angle_results_optimise(save_path): # Create a new .txt file for the results. save_path = os.path.join(save_path, sample.name) with open(os.path.join(save_path, 'optimised_angles.txt'), 'w') as file: - optimiser = Optimiser(sample) # Optimiser for the experiment. + optimiser = Optimiser(sample) # Optimiser for the experiment. # Optimise the experiment using 1-4 angles. for i, num_angles in enumerate([1, 2, 3, 4]): @@ -89,18 +79,24 @@ def _angle_results_optimise(save_path): # Convert to percentages. angles, splits, val = results - splits = np.array(splits)*100 + splits = np.array(splits) * 100 # Round the optimisation function value to 4 significant figures. - val = np.format_float_positional(val, precision=4, unique=False, - fractional=False, trim='k') + val = np.format_float_positional( + val, precision=4, unique=False, fractional=False, trim='k' + ) # Save the conditions, objective value and computation time. - file.write('----------- {} Angles -----------\n'.format(num_angles)) + file.write( + '----------- {} Angles -----------\n'.format(num_angles) + ) file.write('Angles: {}\n'.format(list(np.round(angles, 2)))) file.write('Splits (%): {}\n'.format(list(np.round(splits, 1)))) file.write('Objective value: {}\n'.format(val)) - file.write('Computation time: {}\n\n'.format(round(end-start, 1))) + file.write( + 'Computation time: {}\n\n'.format(round(end - start, 1)) + ) + if __name__ == '__main__': save_path = './results' diff --git a/hogben/contrasts.py b/hogben/contrasts.py index 520143d3..77a5f442 100644 --- a/hogben/contrasts.py +++ b/hogben/contrasts.py @@ -1,5 +1,4 @@ import os -import sys import time import matplotlib.pyplot as plt @@ -17,7 +16,7 @@ def _contrast_results_visualise(save_path): save_path (str): path to directory to save results to. """ - from models.bilayers import BilayerDMPC, BilayerDPPC + from models.bilayers import BilayerDMPC # Choose sample here. bilayer = BilayerDMPC() @@ -25,12 +24,17 @@ def _contrast_results_visualise(save_path): # Number of points and counting times for each angle to simulate. angle_times = [(0.7, 100, 10), (2.3, 100, 40)] - # Visualise single contrast choices assuming different initial measurements. + # Visualise single contrast choices assuming different initial + # measurements. contrast_range = np.linspace(-0.56, 6.36, 500) - contrast_choice_single(bilayer, contrast_range, [], angle_times, save_path, 'initial') - contrast_choice_single(bilayer, contrast_range, [6.36], angle_times, save_path, 'D2O') - contrast_choice_single(bilayer, contrast_range, [-0.56], angle_times, save_path, 'H2O') - contrast_choice_single(bilayer, contrast_range, [-0.56, 6.36], angle_times, save_path, 'D2O_H2O') + contrast_choice_single(bilayer, contrast_range, [], angle_times, + save_path, 'initial') + contrast_choice_single(bilayer, contrast_range, [6.36], angle_times, + save_path, 'D2O') + contrast_choice_single(bilayer, contrast_range, [-0.56], angle_times, + save_path, 'H2O') + contrast_choice_single(bilayer, contrast_range, [-0.56, 6.36], angle_times, + save_path, 'D2O_H2O') # Investigate contrast pair choices assuming no prior measurement. contrast_range = np.linspace(-0.56, 6.36, 75) @@ -40,6 +44,7 @@ def _contrast_results_visualise(save_path): bilayer.nested_sampling([6.36, 6.36], angle_times, save_path, 'D2O_D2O') bilayer.nested_sampling([-0.56, 6.36], angle_times, save_path, 'H2O_D2O') + def _contrast_results_optimise(save_path): """Optimises the choice of contrasts for a bilayer sample. @@ -47,13 +52,13 @@ def _contrast_results_optimise(save_path): save_path (str): path to directory to save results to. """ - from bilayers import BilayerDMPC, BilayerDPPC + from bilayers import BilayerDMPC # Choose sample here. bilayer = BilayerDMPC() # Time budget for experiment. - total_time = 1000 # A large time improves DE convergence. + total_time = 1000 # A large time improves DE convergence. # Points and proportion of total counting time for each angle. angle_splits = [(0.7, 100, 0.2), (2.3, 100, 0.8)] @@ -64,7 +69,7 @@ def _contrast_results_optimise(save_path): # Create a new .txt file for the results. save_path = os.path.join(save_path, bilayer.name) with open(os.path.join(save_path, 'optimised_contrasts.txt'), 'w') as file: - optimiser = Optimiser(bilayer) # Optimiser for the experiment. + optimiser = Optimiser(bilayer) # Optimiser for the experiment. # Optimise the experiment using 1-4 contrasts. for i, num_contrasts in enumerate([1, 2, 3, 4]): @@ -79,18 +84,23 @@ def _contrast_results_optimise(save_path): end = time.time() contrasts, splits, val = results - splits = np.array(splits)*100 # Convert to percentages. + splits = np.array(splits) * 100 # Convert to percentages. # Round the optimisation function value to 4 significant figures. val = np.format_float_positional(val, precision=4, unique=False, fractional=False, trim='k') # Write the conditions, objective value and computation time. - file.write('----------- {} Contrasts -----------\n'.format(num_contrasts)) + file.write( + '----------- {} Contrasts -----------\n'.format(num_contrasts)) file.write('Contrasts: {}\n'.format(list(np.round(contrasts, 2)))) file.write('Splits (%): {}\n'.format(list(np.round(splits, 1)))) file.write('Objective value: {}\n'.format(val)) - file.write('Computation time: {}\n\n'.format(round(end-start, 1))) + file.write( + 'Computation time: {}\n\n'.format( + round( + end - start, 1))) + def _figure_2(save_path): """Creates figure 2 for the paper. The figure shows the minimum eigenvalue @@ -117,36 +127,36 @@ def _figure_2(save_path): h2o_split_2 = 0.245 # Define how much the third contrast will be measured for. - nxt_split_1 = 1-d2o_split_1-h2o_split_1 - nxt_split_2 = 1-d2o_split_2-h2o_split_2 + nxt_split_1 = 1 - d2o_split_1 - h2o_split_1 + nxt_split_2 = 1 - d2o_split_2 - h2o_split_2 # Counting times for each angle when measuring D2O. - d2o_angle_times_1 = [(angle, points, total_time*split*d2o_split_1) + d2o_angle_times_1 = [(angle, points, total_time * split * d2o_split_1) for angle, points, split in angle_splits] - d2o_angle_times_2 = [(angle, points, total_time*split*d2o_split_2) + d2o_angle_times_2 = [(angle, points, total_time * split * d2o_split_2) for angle, points, split in angle_splits] # Counting times for each angle when measuring H2O. - h2o_angle_times_1 = [(angle, points, total_time*split*h2o_split_1) + h2o_angle_times_1 = [(angle, points, total_time * split * h2o_split_1) for angle, points, split in angle_splits] - h2o_angle_times_2 = [(angle, points, total_time*split*h2o_split_2) + h2o_angle_times_2 = [(angle, points, total_time * split * h2o_split_2) for angle, points, split in angle_splits] # Counting times for each angle when measuring the third contrast. - nxt_angle_times_1 = [(angle, points, total_time*split*(nxt_split_1)) + nxt_angle_times_1 = [(angle, points, total_time * split * (nxt_split_1)) for angle, points, split in angle_splits] - nxt_angle_times_2 = [(angle, points, total_time*split*(nxt_split_2)) + nxt_angle_times_2 = [(angle, points, total_time * split * (nxt_split_2)) for angle, points, split in angle_splits] # Information from the D2O and H2O contrasts for each model. - g_init_1 = (sample_1.angle_info(d2o_angle_times_1, [6.36]) + - sample_1.angle_info(h2o_angle_times_1, [-0.56])) + g_init_1 = (sample_1.angle_info(d2o_angle_times_1, [6.36]) + + sample_1.angle_info(h2o_angle_times_1, [-0.56])) - g_init_2 = (sample_2.angle_info(d2o_angle_times_2, [6.36]) + - sample_2.angle_info(h2o_angle_times_2, [-0.56])) + g_init_2 = (sample_2.angle_info(d2o_angle_times_2, [6.36]) + + sample_2.angle_info(h2o_angle_times_2, [-0.56])) # Calculate the minimum eigenvalue for each choice of third contrast. min_eigs_1, min_eigs_2 = [], [] @@ -156,8 +166,8 @@ def _figure_2(save_path): g_new_1 = sample_1.contrast_info(nxt_angle_times_1, [new_contrast]) g_new_2 = sample_2.contrast_info(nxt_angle_times_2, [new_contrast]) - min_eigs_1.append(np.linalg.eigvalsh(g_init_1+g_new_1)[0]) - min_eigs_2.append(np.linalg.eigvalsh(g_init_2+g_new_2)[0]) + min_eigs_1.append(np.linalg.eigvalsh(g_init_1 + g_new_1)[0]) + min_eigs_2.append(np.linalg.eigvalsh(g_init_2 + g_new_2)[0]) # Create the plot of third contrast versus minimum eigenvalue. fig = plt.figure() @@ -174,11 +184,12 @@ def _figure_2(save_path): ax1.set_xlabel(x_label, fontsize=11, weight='bold') ax1.set_ylabel(y_label, fontsize=11, weight='bold', color='b') ax2.set_ylabel(y_label, fontsize=11, weight='bold', color='g') - ax1.legend(line_1+line_2, ['DMPC', 'DPPC/RaLPS'], loc=0) + ax1.legend(line_1 + line_2, ['DMPC', 'DPPC/RaLPS'], loc=0) # Save the plot. save_plot(fig, save_path, 'figure_2') + if __name__ == '__main__': save_path = './results' _contrast_results_visualise(save_path) diff --git a/hogben/data/directbeams/make_beam_spectra.py b/hogben/data/directbeams/make_beam_spectra.py index da264e68..61f8945d 100644 --- a/hogben/data/directbeams/make_beam_spectra.py +++ b/hogben/data/directbeams/make_beam_spectra.py @@ -1,86 +1,217 @@ -#Load files and do flood correction +# Load files and do flood correction # import mantid algorithms, numpy and matplotlib from mantid.simpleapi import * -import matplotlib.pyplot as plt -import numpy as np -def pre_process(runno,detectors='75-95'): - LoadISISNexus(str(runno),OutputWorkspace=str(runno),LoadMonitors="Exclude") - ConvertUnits(InputWorkspace=str(runno), OutputWorkspace=str(runno), Target='Wavelength', AlignBins=True) + +def pre_process(runno, detectors='75-95'): + LoadISISNexus( + str(runno), OutputWorkspace=str(runno), LoadMonitors='Exclude' + ) + ConvertUnits( + InputWorkspace=str(runno), + OutputWorkspace=str(runno), + Target='Wavelength', + AlignBins=True, + ) try: - SumSpectra(InputWorkspace=str(runno), OutputWorkspace=str(runno), ListOfWorkspaceIndices=detectors) - except: - ExtractSpectra(InputWorkspace=str(runno), OutputWorkspace=str(runno), WorkspaceIndexList=detectors) - ReplaceSpecialValues(InputWorkspace=str(runno), OutputWorkspace=str(runno), NaNValue=0, InfinityValue=0) - NormaliseByCurrent(str(runno), OutputWorkspace = str(runno)) - CropWorkspace(str(runno),XMin=1.0, XMax=17.0, OutputWorkspace = str(runno)) - Rebin(str(runno), Params = 0.005, OutputWorkspace = str(runno)) + SumSpectra( + InputWorkspace=str(runno), + OutputWorkspace=str(runno), + ListOfWorkspaceIndices=detectors, + ) + except BaseException: + ExtractSpectra( + InputWorkspace=str(runno), + OutputWorkspace=str(runno), + WorkspaceIndexList=detectors, + ) + ReplaceSpecialValues( + InputWorkspace=str(runno), + OutputWorkspace=str(runno), + NaNValue=0, + InfinityValue=0, + ) + NormaliseByCurrent(str(runno), OutputWorkspace=str(runno)) + CropWorkspace(str(runno), XMin=1.0, XMax=17.0, OutputWorkspace=str(runno)) + Rebin(str(runno), Params=0.005, OutputWorkspace=str(runno)) -#INTER -pre_process(62096) # Both DB's to get correct flux 0.5 degrees -pre_process(62097) -Scale("62096",Factor=13,OutputWorkspace="62096") #Quickest way is to do this by hand I think -a=0.7/0.5 -Scale("62096",Factor=a,OutputWorkspace="62096") #To make it effectively at 0.7 degrees, same as OFFSPEC beam -SaveAscii("62096",Filename="/mnt/ceph/home/jc15575/Desktop/Enough/enough-is-enough-main/simulation/INTER_Air_DB.dat",WriteSpectrumID=False,ColumnHeader=False) -#POLREF -pre_process(34365,detectors='272-285') #0.5 degrees I think -a=0.7/0.5 -Scale("34365",Factor=a,OutputWorkspace="34365") #To make it effectively at 0.7 degrees, same as OFFSPEC beam -SaveAscii("34365",Filename="/mnt/ceph/home/jc15575/Desktop/Enough/enough-is-enough-main/simulation/POLREF_Air_DB.dat",WriteSpectrumID=False,ColumnHeader=False) +# INTER +pre_process(62096) # Both DB's to get correct flux 0.5 degrees +pre_process(62097) +Scale( + '62096', Factor=13, OutputWorkspace='62096' +) # Quickest way is to do this by hand I think +a = 0.7 / 0.5 +Scale( + '62096', Factor=a, OutputWorkspace='62096' +) # To make it effectively at 0.7 degrees, same as OFFSPEC beam +SaveAscii( + '62096', + Filename=('/mnt/ceph/home/jc15575/Desktop/Enough/enough-is-enough-main' + '/simulation/INTER_Air_DB.dat'), + WriteSpectrumID=False, + ColumnHeader=False, +) -#OFFSPEC -pre_process(62802,detectors='380-420') # Both DB's to get correct flux, 0.7 degrees -pre_process(62803,detectors='380-420') -Scale("62802",Factor=3.7,OutputWorkspace="62802") #Quickest way is to do this by hand I think -SaveAscii("62802",Filename="/mnt/ceph/home/jc15575/Desktop/Enough/enough-is-enough-main/simulation/OFFSPEC_Air_DB.dat",WriteSpectrumID=False,ColumnHeader=False) -pre_process(53439,detectors='380-420') #Polarised +# POLREF +pre_process(34365, detectors='272-285') # 0.5 degrees I think +a = 0.7 / 0.5 +Scale( + '34365', Factor=a, OutputWorkspace='34365' +) # To make it effectively at 0.7 degrees, same as OFFSPEC beam +SaveAscii( + '34365', + Filename=('/mnt/ceph/home/jc15575/Desktop/Enough/enough-is-enough-main/' + 'simulation/POLREF_Air_DB.dat'), + WriteSpectrumID=False, + ColumnHeader=False, +) -#SURF -pre_process(137344,detectors=0) #0.35 degrees I think -a=0.7/0.35 -Scale("137344",Factor=a,OutputWorkspace="137344") -SaveAscii("137344",Filename="/mnt/ceph/home/jc15575/Desktop/Enough/enough-is-enough-main/simulation/SURF_Air_DB.dat",WriteSpectrumID=False,ColumnHeader=False) +# OFFSPEC +pre_process( + 62802, detectors='380-420' +) # Both DB's to get correct flux, 0.7 degrees +pre_process(62803, detectors='380-420') +Scale( + '62802', Factor=3.7, OutputWorkspace='62802' +) # Quickest way is to do this by hand I think +SaveAscii( + '62802', + Filename=('/mnt/ceph/home/jc15575/Desktop/Enough/enough-is-enough-main' + '/simulation/OFFSPEC_Air_DB.dat'), + WriteSpectrumID=False, + ColumnHeader=False, +) +pre_process(53439, detectors='380-420') # Polarised +# SURF +pre_process(137344, detectors=0) # 0.35 degrees I think +a = 0.7 / 0.35 +Scale('137344', Factor=a, OutputWorkspace='137344') +SaveAscii( + '137344', + Filename=('/mnt/ceph/home/jc15575/Desktop/Enough/enough-is-enough-main' + '/simulation/SURF_Air_DB.dat'), + WriteSpectrumID=False, + ColumnHeader=False, +) -LoadNexusProcessed(Filename=r'C:/Users/npi34092/Dropbox/Offspec_reduction\GC_Flood_2-19.nxs', OutputWorkspace='wsf1b_flood') -ConvertUnits(InputWorkspace='wsf1b_flood', OutputWorkspace='wsf1b_flood', Target='Wavelength', AlignBins=True) -LoadISISNexus(Filename=r'\\isis.cclrc.ac.uk\inst$\ndxoffspec\instrument\data\cycle_18_4\OFFSPEC00050612.nxs', OutputWorkspace='wtemp') +LoadNexusProcessed( + Filename=r'C:/Users/npi34092/Dropbox/Offspec_reduction\GC_Flood_2-19.nxs', + OutputWorkspace='wsf1b_flood', +) +ConvertUnits( + InputWorkspace='wsf1b_flood', + OutputWorkspace='wsf1b_flood', + Target='Wavelength', + AlignBins=True, +) +LoadISISNexus( + Filename=(r'\\isis.cclrc.ac.uk\inst$\ndxoffspec\instrument\data' + r'\cycle_18_4\OFFSPEC00050612.nxs'), + OutputWorkspace='wtemp', +) CloneWorkspace(InputWorkspace='wtemp', OutputWorkspace='w55') -LoadISISNexus(Filename=r'\\isis.cclrc.ac.uk\inst$\ndxoffspec\instrument\data\cycle_18_4\OFFSPEC00050613.nxs', OutputWorkspace='wtemp') +LoadISISNexus( + Filename=(r'\\isis.cclrc.ac.uk\inst$\ndxoffspec\instrument\data' + r'\cycle_18_4\OFFSPEC00050613.nxs'), + OutputWorkspace='wtemp', +) CloneWorkspace(InputWorkspace='wtemp', OutputWorkspace='w56') -ConvertUnits(InputWorkspace='w55', OutputWorkspace='w55', Target='Wavelength', AlignBins=True) -CropWorkspace(InputWorkspace='w55', OutputWorkspace='w55det', StartWorkspaceIndex=5, EndWorkspaceIndex=771) -RebinToWorkspace(WorkspaceToRebin='wsf1b_flood', WorkspaceToMatch='w55det', OutputWorkspace='wsf1b_flood_reb') -Divide(LHSWorkspace='w55det', RHSWorkspace='wsf1b_flood_reb', OutputWorkspace='w55det', AllowDifferentNumberSpectra=True) +ConvertUnits( + InputWorkspace='w55', + OutputWorkspace='w55', + Target='Wavelength', + AlignBins=True, +) +CropWorkspace( + InputWorkspace='w55', + OutputWorkspace='w55det', + StartWorkspaceIndex=5, + EndWorkspaceIndex=771, +) +RebinToWorkspace( + WorkspaceToRebin='wsf1b_flood', + WorkspaceToMatch='w55det', + OutputWorkspace='wsf1b_flood_reb', +) +Divide( + LHSWorkspace='w55det', + RHSWorkspace='wsf1b_flood_reb', + OutputWorkspace='w55det', + AllowDifferentNumberSpectra=True, +) -#Sum up the direct beam spectra and do the same for the second workspace -SumSpectra(InputWorkspace='w55det', OutputWorkspace='w55norm', ListOfWorkspaceIndices='384-414') -ReplaceSpecialValues(InputWorkspace='w55norm', OutputWorkspace='w55norm', NaNValue=0, InfinityValue=0) -ConvertUnits(InputWorkspace='w56', OutputWorkspace='w56', Target='Wavelength', AlignBins=True) -CropWorkspace(InputWorkspace='w56', OutputWorkspace='w56det', StartWorkspaceIndex=5, EndWorkspaceIndex=771) -RebinToWorkspace(WorkspaceToRebin='wsf1b_flood', WorkspaceToMatch='w56det', OutputWorkspace='wsf1b_flood_reb') -Divide(LHSWorkspace='w56det', RHSWorkspace='wsf1b_flood_reb', OutputWorkspace='w56det', AllowDifferentNumberSpectra=True) -SumSpectra(InputWorkspace='w56det', OutputWorkspace='w56norm', ListOfWorkspaceIndices='384-414') -ReplaceSpecialValues(InputWorkspace='w56norm', OutputWorkspace='w56norm', NaNValue=0, InfinityValue=0) +# Sum up the direct beam spectra and do the same for the second workspace +SumSpectra( + InputWorkspace='w55det', + OutputWorkspace='w55norm', + ListOfWorkspaceIndices='384-414', +) +ReplaceSpecialValues( + InputWorkspace='w55norm', + OutputWorkspace='w55norm', + NaNValue=0, + InfinityValue=0, +) +ConvertUnits( + InputWorkspace='w56', + OutputWorkspace='w56', + Target='Wavelength', + AlignBins=True, +) +CropWorkspace( + InputWorkspace='w56', + OutputWorkspace='w56det', + StartWorkspaceIndex=5, + EndWorkspaceIndex=771, +) +RebinToWorkspace( + WorkspaceToRebin='wsf1b_flood', + WorkspaceToMatch='w56det', + OutputWorkspace='wsf1b_flood_reb', +) +Divide( + LHSWorkspace='w56det', + RHSWorkspace='wsf1b_flood_reb', + OutputWorkspace='w56det', + AllowDifferentNumberSpectra=True, +) +SumSpectra( + InputWorkspace='w56det', + OutputWorkspace='w56norm', + ListOfWorkspaceIndices='384-414', +) +ReplaceSpecialValues( + InputWorkspace='w56norm', + OutputWorkspace='w56norm', + NaNValue=0, + InfinityValue=0, +) -#For a single direct beam I think this will do. s1/s2(h) = 1.95 (30) / 0.1725 (30) -NormaliseByCurrent('w56norm', OutputWorkspace = 'w56norm') -CropWorkspace('w56norm',XMin=1.0, XMax=14.0, OutputWorkspace = 'w56norm') -Rebin('w56norm', Params = 0.05, OutputWorkspace = 'w56norm') -''' +# For a single direct beam I think this will do. s1/s2(h) = 1.95 (30) / +# 0.1725 (30) +NormaliseByCurrent('w56norm', OutputWorkspace='w56norm') +CropWorkspace('w56norm', XMin=1.0, XMax=14.0, OutputWorkspace='w56norm') +Rebin('w56norm', Params=0.05, OutputWorkspace='w56norm') +""" #Scale the direct beam with smaller slits correctly and add them together -Integration(InputWorkspace='w55norm', OutputWorkspace='w55int', RangeLower=1, RangeUpper=14) -Integration(InputWorkspace='w56norm', OutputWorkspace='w56int', RangeLower=1, RangeUpper=14) -Multiply(LHSWorkspace='w55norm', RHSWorkspace='w56int', OutputWorkspace='w55norm') -Divide(LHSWorkspace='w55norm', RHSWorkspace='w55int', OutputWorkspace='w55norm') +Integration(InputWorkspace='w55norm', OutputWorkspace='w55int', RangeLower=1, + RangeUpper=14) +Integration(InputWorkspace='w56norm', OutputWorkspace='w56int', RangeLower=1, + RangeUpper=14) +Multiply(LHSWorkspace='w55norm', RHSWorkspace='w56int', + OutputWorkspace='w55norm') +Divide(LHSWorkspace='w55norm', RHSWorkspace='w55int', + OutputWorkspace='w55norm') MultiplyRange(InputWorkspace='w56norm', OutputWorkspace='w56norm', EndBin=157) -WeightedMean(InputWorkspace1='w55norm', InputWorkspace2='w56norm', OutputWorkspace='DBair03') +WeightedMean(InputWorkspace1='w55norm', InputWorkspace2='w56norm', + OutputWorkspace='DBair03') # import mantid algorithms from mantid.simpleapi import * -''' \ No newline at end of file +""" diff --git a/hogben/kinetics.py b/hogben/kinetics.py index 7c60cecb..3e4f8e22 100644 --- a/hogben/kinetics.py +++ b/hogben/kinetics.py @@ -1,6 +1,5 @@ import numpy as np import os -import sys import matplotlib.pyplot as plt @@ -36,18 +35,18 @@ def _kinetics_results_visualise(save_path): # Iterate over each contrast and angle being considered. x, y, infos = [], [], [] - n = len(angle_range)*len(contrast_range) # Number of calculations. + n = len(angle_range) * len(contrast_range) # Number of calculations. for i, contrast_sld in enumerate(contrast_range): # Display progress. if i % 5 == 0: - print('>>> {0}/{1}'.format(i*len(angle_range), n)) + print('>>> {0}/{1}'.format(i * len(angle_range), n)) for angle in angle_range: # Record the "true" lipid APM value. apm = monolayer.lipid_apm.value # Split the time budget based on number of APM values. - angle_times = [(angle, points, time/len(apm_range))] + angle_times = [(angle, points, time / len(apm_range))] # Calculate the lipid APM Fisher information for each value. information = 0 @@ -57,13 +56,13 @@ def _kinetics_results_visualise(save_path): g = monolayer.contrast_info(angle_times, [contrast_sld]) information += g[0, 0] - monolayer.lipid_apm.value = apm # Reset the APM parameter. + monolayer.lipid_apm.value = apm # Reset the APM parameter. infos.append(information) x.append(contrast_sld) y.append(angle) # Create the plot of angle and contrast SLD versus Fisher information. - fig = plt.figure(figsize=[10,8]) + fig = plt.figure(figsize=[10, 8]) ax = fig.add_subplot(111, projection='3d') # Create the surface plot and add colour bar. @@ -77,7 +76,7 @@ def _kinetics_results_visualise(save_path): ax.set_xlabel(x_label, fontsize=11, weight='bold') ax.set_ylabel(y_label, fontsize=11, weight='bold') ax.set_zlabel(z_label, fontsize=11, weight='bold') - ax.ticklabel_format(axis='z', style='sci', scilimits=(0,0)) + ax.ticklabel_format(axis='z', style='sci', scilimits=(0, 0)) # Save the plot. save_path = os.path.join(save_path, monolayer.name) @@ -100,6 +99,7 @@ def _kinetics_results_visualise(save_path): maximum = np.argmax(infos) return x[maximum], y[maximum] + def _kinetics_results_optimise(save_path): """Optimises the choice of measurement angle and contrast SLD for the DPPG monolayer model where the lipid area per molecule is decreasing @@ -110,7 +110,7 @@ def _kinetics_results_optimise(save_path): """ # Counting time and number of points to use when simulating. - time = 1000 # A large time improves DE convergence. + time = 1000 # A large time improves DE convergence. points = 100 # Intervals containing angles and contrasts to consider. @@ -133,7 +133,7 @@ def _kinetics_results_optimise(save_path): # Optimise angle and contrast. res = differential_evolution(_optimisation_func, bounds, - args=[monolayer]+args, polish=False, + args=[monolayer] + args, polish=False, tol=0.001, updating='deferred', workers=-1, disp=False) @@ -145,6 +145,7 @@ def _kinetics_results_optimise(save_path): file.write('Contrast: {}\n'.format(round(contrast, 2))) file.write('Objective value: {}\n\n'.format(res.fun)) + def _optimisation_func(x, sample, apm_range, points, time): """Defines the function for optimising the DPPG monolayer measurement angle and contrast. @@ -163,7 +164,7 @@ def _optimisation_func(x, sample, apm_range, points, time): """ # Define the points and counting times for each angle to simulate. angle, contrast_sld = x[0], x[1] - angle_times = [(angle, points, time/len(apm_range))] + angle_times = [(angle, points, time / len(apm_range))] # Calculate the information in the experiment using the given conditions. information = 0 @@ -176,6 +177,7 @@ def _optimisation_func(x, sample, apm_range, points, time): sample.lipid_apm.value = apm return -information + if __name__ == '__main__': save_path = './results' diff --git a/hogben/magnetism.py b/hogben/magnetism.py index 20e19b7c..b99c67cd 100644 --- a/hogben/magnetism.py +++ b/hogben/magnetism.py @@ -1,6 +1,5 @@ import matplotlib.pyplot as plt import os -import sys import numpy as np @@ -26,33 +25,32 @@ def _magnetism_results_visualise(save_path): sample.pt_mag.value = 0.01638 # Number of points and counting times for each angle to simulate. - angle_times = [(0.5, 100, 20), - (1.0, 100, 40), - (2.0, 100, 80)] + angle_times = [(0.5, 100, 20), (1.0, 100, 40), (2.0, 100, 80)] # Range of YIG and Pt layer thicknesses to calculate over. yig_thick_range = np.linspace(400, 900, 75) - pt_thick_range = np.concatenate((np.linspace(21.5, 30, 50), - np.linspace(30, 100, 50))) + pt_thick_range = np.concatenate( + (np.linspace(21.5, 30, 50), np.linspace(30, 100, 50)) + ) # Iterate over each YIG and Pt layer thickness being considered. x, y, infos = [], [], [] - n = len(pt_thick_range)*len(yig_thick_range) # Number of calculations. + n = len(pt_thick_range) * len(yig_thick_range) # Number of calculations. for i, yig_thick in enumerate(yig_thick_range): # Display progress. if i % 5 == 0: - print('>>> {0}/{1}'.format(i*len(pt_thick_range), n)) + print('>>> {0}/{1}'.format(i * len(pt_thick_range), n)) for pt_thick in pt_thick_range: # Calculate the Fisher information using current thicknesses. g = sample.underlayer_info(angle_times, yig_thick, pt_thick) - infos.append(g[0,0]) + infos.append(g[0, 0]) x.append(yig_thick) y.append(pt_thick) # Create plot of YIG and Pt layer thicknesses versus Fisher information. - fig = plt.figure(figsize=[10,8]) + fig = plt.figure(figsize=[10, 8]) ax = fig.add_subplot(111, projection='3d') # Create the surface plot and add colour bar. @@ -66,7 +64,7 @@ def _magnetism_results_visualise(save_path): ax.set_xlabel(x_label, fontsize=11, weight='bold') ax.set_ylabel(y_label, fontsize=11, weight='bold') ax.set_zlabel(z_label, fontsize=11, weight='bold') - ax.ticklabel_format(axis='z', style='sci', scilimits=(0,0)) + ax.ticklabel_format(axis='z', style='sci', scilimits=(0, 0)) # Save the plot. save_path = os.path.join(save_path, sample.name) @@ -85,6 +83,7 @@ def _magnetism_results_visualise(save_path): maximum = np.argmax(infos) return x[maximum], y[maximum] + def _magnetism_results_optimise(save_path): """Optimises the choice of YIG and Pt layer thicknesses for the magnetic YIG sample. @@ -97,9 +96,7 @@ def _magnetism_results_optimise(save_path): sample.pt_mag.value = 0.01638 # Number of points and counting times for each angle to simulate. - angle_times = [(0.5, 100, 20), - (1.0, 100, 40), - (2.0, 100, 80)] + angle_times = [(0.5, 100, 20), (1.0, 100, 40), (2.0, 100, 80)] # Itervals of YIG and Pt layer thicknesses to optimise over. yig_thick_bounds = (200, 900) @@ -114,10 +111,16 @@ def _magnetism_results_optimise(save_path): file_path = os.path.join(save_path, 'optimised_underlayers.txt') with open(file_path, 'w') as file: # Optimise thicknesses using differential evolution. - res = differential_evolution(_optimisation_func, bounds, args=args, - polish=False, tol=0.001, - updating='deferred', workers=-1, - disp=False) + res = differential_evolution( + _optimisation_func, + bounds, + args=args, + polish=False, + tol=0.001, + updating='deferred', + workers=-1, + disp=False, + ) yig_thick, pt_thick = res.x[0], res.x[1] @@ -132,6 +135,7 @@ def _magnetism_results_optimise(save_path): file.write('Pt Thickness: {}\n'.format(round(pt_thick, 1))) file.write('Fisher Information: {}\n\n'.format(-res.fun)) + def _optimisation_func(x, sample, angle_times): """Defines the function for optimising YIG and Pt layer thicknesses for the magnetic YIG sample. @@ -148,7 +152,8 @@ def _optimisation_func(x, sample, angle_times): """ # Calculate the Fisher information matrix using the given thicknesses. g = sample.underlayer_info(angle_times, x[0], x[1]) - return -g[0,0] + return -g[0, 0] + def _magnetism_results_ratios(save_path): """Calculates log ratio of likelihoods between two models, one with an @@ -167,28 +172,28 @@ def _magnetism_results_ratios(save_path): times = np.linspace(40, 4000, 250) # Points and split of total counting time for each angle to simulate. - angle_splits = [(0.5, 100, 1/7), - (1.0, 100, 2/7), - (2.0, 100, 4/7)] + angle_splits = [(0.5, 100, 1 / 7), (1.0, 100, 2 / 7), (2.0, 100, 4 / 7)] # Calculate log ratio of likelihoods with optimal and sub-optimal designs. _calc_log_ratios(26, times, angle_splits, save_path) _calc_log_ratios(80, times, angle_splits, save_path) # Create the plot of counting time versus log ratio of likelihoods. - fig = plt.figure(figsize=(6,7)) + fig = plt.figure(figsize=(6, 7)) ax = fig.add_subplot(111) # Iterate over the results for the two designs. save_path = os.path.join(save_path, 'YIG_sample') for pt_thick in [26, 80]: - file_path = os.path.join(save_path, 'log_ratios_{}.csv'.format(pt_thick)) + file_path = os.path.join( + save_path, 'log_ratios_{}.csv'.format(pt_thick) + ) # Load and plot the calculated ratios. data = np.loadtxt(file_path, delimiter=',') - times, ratios = data[:,0], data[:,1] + times, ratios = data[:, 0], data[:, 1] - ax.plot(1.5*times, ratios, label='{}Å Pt Thickness'.format(pt_thick)) + ax.plot(1.5 * times, ratios, label='{}Å Pt Thickness'.format(pt_thick)) ax.set_xlabel('Counting Time (min.)', fontsize=11, weight='bold') ax.set_ylabel('Log Ratio of Likelihoods', fontsize=11, weight='bold') @@ -197,6 +202,7 @@ def _magnetism_results_ratios(save_path): # Save the plot. save_plot(fig, save_path, 'log_ratios') + def _calc_log_ratios(pt_thick, times, angle_splits, save_path): """Calculates log ratio of likelihoods between two models, one with an induced moment in the YIG sample Pt layer and one with no moment, as @@ -219,8 +225,10 @@ def _calc_log_ratios(pt_thick, times, angle_splits, save_path): # Get the ratio for 100 simulated data sets using the time. for _ in range(100): # Define the number of points and times for each angle. - angle_times = [(angle, points, split*total_time) - for angle, points, split in angle_splits] + angle_times = [ + (angle, points, split * total_time) + for angle, points, split in angle_splits + ] # Simulate data for the YIG sample with a 1 uB/atom magnetic # moment in the Pt layer. @@ -243,7 +251,7 @@ def _calc_log_ratios(pt_thick, times, angle_splits, save_path): logl_2 = _logl(models) # Record the ratio of likelihoods. - ratio = logl_1-logl_2 + ratio = logl_1 - logl_2 ratios.append(ratio) # Calculate and save median ratio. @@ -251,6 +259,7 @@ def _calc_log_ratios(pt_thick, times, angle_splits, save_path): file.write('{0},{1}\n'.format(total_time, median_ratio)) print(median_ratio) + def _logl(models): """Calculates the log-likelihood for a given list of `models` corresponding to simulated spin states. @@ -263,7 +272,7 @@ def _logl(models): """ # Extract the Q, R, dR and model R for the simulated spin states. - q, r, dr, r_model = [], [] , [], [] + q, r, dr, r_model = [], [], [], [] for model in models: probe = model.probe.xs[model.probe.spin_state] q.append(probe.Q) @@ -278,7 +287,8 @@ def _logl(models): r_model = np.concatenate(r_model) # Calculate the log-likelihood over all the models. - return -0.5*np.sum(((r-r_model)/dr)**2 + np.log(2*np.pi*dr**2)) + return -0.5 * np.sum((r - r_model) / dr) ** 2 + np.log(2 * np.pi * dr**2) + if __name__ == '__main__': save_path = './results' diff --git a/hogben/models/base.py b/hogben/models/base.py index 2b9beb1e..9e8bef07 100644 --- a/hogben/models/base.py +++ b/hogben/models/base.py @@ -1,7 +1,6 @@ import os -import sys - from abc import ABC, abstractmethod +from typing import Optional import matplotlib.pyplot as plt @@ -19,33 +18,41 @@ class VariableAngle(ABC): """Abstract class representing whether the measurement angle of a sample can be varied.""" + @abstractmethod def angle_info(self): + """Calculates the Fisher information matrix for a sample measured - over a number of angles.""" + over a number of angles.""" pass + class VariableContrast(ABC): """Abstract class representing whether the contrast of a sample - can be varied.""" + dan be varied.""" + @abstractmethod def contrast_info(self): """Calculates the Fisher information matrix for a sample with contrasts measured over a number of angles.""" pass + class VariableUnderlayer(ABC): """Abstract class representing whether the underlayer(s) of a sample can be varied.""" + @abstractmethod def underlayer_info(self): """Calculates the Fisher information matrix for a sample with - underlayers, and contrasts measured over a number of angles.""" + underlayers, and contrasts measured over a number of angles.""" pass + class BaseSample(VariableAngle): """Abstract class representing a "standard" neutron reflectometry sample - defined by a series of contiguous layers.""" + defined by a series of contiguous layers.""" + @abstractmethod def sld_profile(self): """Plots the SLD profile of the sample.""" @@ -61,10 +68,12 @@ def nested_sampling(self): """Runs nested sampling on measured or simulated data of the sample.""" pass + class BaseLipid(BaseSample, VariableContrast, VariableUnderlayer): """Abstract class representing the base class for a lipid model.""" + def __init__(self): - self._create_objectives() # Load experimentally-measured data. + self._create_objectives() # Load experimentally-measured data. @abstractmethod def _create_objectives(self): @@ -133,10 +142,14 @@ def __conditions_info(self, angle_times, contrasts, underlayers): # Simulate data for the contrast. sample = self._using_conditions(contrast, underlayers) contrast_point = (contrast + 0.56) / (6.35 + 0.56) - background_level = 2e-6*contrast_point + 4e-6*(1-contrast_point) - model, data = simulate(sample, angle_times, scale=1, bkg=background_level, dq=2) - qs.append(data[:,0]) - counts.append(data[:,3]) + background_level = (2e-6 * contrast_point + + 4e-6 * (1 - contrast_point) + ) + model, data = simulate( + sample, angle_times, scale=1, bkg=background_level, dq=2 + ) + qs.append(data[:, 0]) + counts.append(data[:, 3]) models.append(model) # Exclude certain parameters if underlayers are being used. @@ -150,8 +163,11 @@ def _using_conditions(self): """Creates a structure describing the given measurement conditions.""" pass - def sld_profile(self, save_path, filename='sld_profile', - ylim=None, legend=True): + def sld_profile(self, + save_path: str, + filename: str = 'sld_profile', + ylim: Optional[tuple] = None, + legend: bool = True) -> None: """Plots the SLD profile of the lipid sample. Args: @@ -186,7 +202,9 @@ def sld_profile(self, save_path, filename='sld_profile', save_path = os.path.join(save_path, self.name) save_plot(fig, save_path, filename) - def reflectivity_profile(self, save_path, filename='reflectivity_profile'): + def reflectivity_profile(self, + save_path: str, + filename: str = 'reflectivity_profile') -> None: """Plots the reflectivity profile of the lipid sample. Args: @@ -205,7 +223,7 @@ def reflectivity_profile(self, save_path, filename='reflectivity_profile'): r_model = objective.model(q) # Offset the data, for clarity. - offset = 10**(-2*i) + offset = 10 ** (-2 * i) r *= offset dr *= offset r_model *= offset @@ -213,7 +231,7 @@ def reflectivity_profile(self, save_path, filename='reflectivity_profile'): # Add the offset in the label. label = self.labels[i] if offset != 1: - label += ' $\mathregular{(x10^{-'+str(2*i)+'})}$' + label += ' $\\mathregular{(x10^{-' + str(2 * i) + '})}$' # Plot the measured data and the model reflectivity. ax.errorbar(q, r, dr, @@ -221,7 +239,7 @@ def reflectivity_profile(self, save_path, filename='reflectivity_profile'): color=colours[i], label=label) ax.plot(q, r_model, color=colours[i], zorder=20) - x_label = '$\mathregular{Q\ (Å^{-1})}$' + x_label = '$\\mathregular{Q\\ (Å^{-1})}$' y_label = 'Reflectivity (arb.)' ax.set_xlabel(x_label, fontsize=11, weight='bold') @@ -235,8 +253,13 @@ def reflectivity_profile(self, save_path, filename='reflectivity_profile'): save_path = os.path.join(save_path, self.name) save_plot(fig, save_path, filename) - def nested_sampling(self, contrasts, angle_times, save_path, filename, - underlayers=None, dynamic=False): + def nested_sampling(self, + contrasts: list, + angle_times: list, + save_path: str, + filename: str, + underlayers=None, + dynamic=False) -> None: """Runs nested sampling on simulated data of the lipid sample. Args: @@ -254,9 +277,15 @@ def nested_sampling(self, contrasts, angle_times, save_path, filename, # Simulate an experiment using the given contrast. sample = self._using_conditions(contrast, underlayers) contrast_point = (contrast + 0.56) / (6.35 + 0.56) - background_level = 2e-6*contrast_point + 4e-6*(1-contrast_point) - model, data = simulate(sample, angle_times, scale=1, bkg=background_level, dq=2) - dataset = refnx.dataset.ReflectDataset([data[:,0], data[:,1], data[:,2]]) + background_level = 2e-6 * contrast_point + 4e-6 * ( + 1 - contrast_point + ) + model, data = simulate( + sample, angle_times, scale=1, bkg=background_level, dq=2 + ) + dataset = refnx.dataset.ReflectDataset( + [data[:, 0], data[:, 1], data[:, 2]] + ) objectives.append(refnx.analysis.Objective(model, dataset)) # Combine objectives into a single global objective. @@ -266,7 +295,9 @@ def nested_sampling(self, contrasts, angle_times, save_path, filename, if underlayers is None: global_objective.varying_parameters = lambda: self.params else: - global_objective.varying_parameters = lambda: self.underlayer_params + global_objective.varying_parameters = ( + lambda: self.underlayer_params + ) # Sample the objective using nested sampling. sampler = Sampler(global_objective) @@ -274,4 +305,4 @@ def nested_sampling(self, contrasts, angle_times, save_path, filename, # Save the sampling corner plot. save_path = os.path.join(save_path, self.name) - save_plot(fig, save_path, 'nested_sampling_'+filename) + save_plot(fig, save_path, 'nested_sampling_' + filename) diff --git a/hogben/models/bilayers.py b/hogben/models/bilayers.py index e0030d4a..9eaae1a9 100644 --- a/hogben/models/bilayers.py +++ b/hogben/models/bilayers.py @@ -22,12 +22,13 @@ def neutron_scattering_length(formula: str): Determine the neutron scattering length for a chemical formula. :param formula: Chemical formula. - :return: Real and imaginary descriptors for the scattering length in angstrom. + :return: Real and imaginary descriptors for the scattering length in + Ångström. """ formula_as_dict = parse_formula(formula) scattering_length = 0 + 0j for key, value in formula_as_dict.items(): - scattering_length += (pt.elements.symbol(key).neutron.b_c * value) + scattering_length += pt.elements.symbol(key).neutron.b_c * value if pt.elements.symbol(key).neutron.b_c_i: inc = pt.elements.symbol(key).neutron.b_c_i else: @@ -69,12 +70,12 @@ class BilayerPOPC(BaseLipid): objectives (list): objectives corresponding to each measured contrast. """ + def __init__(self): self.name = 'POPC_bilayer' - self.data_path = os.path.join(os.path.dirname(__file__), - '..', - 'data', - 'POPC_bilayer') + self.data_path = os.path.join( + os.path.dirname(__file__), '..', 'data', 'POPC_bilayer' + ) self.labels = ['Si-D2O', 'Si-POPC-D2O', 'Si-POPC-H2O'] self.distances = np.linspace(-20, 95, 500) @@ -83,46 +84,66 @@ def __init__(self): self.dq = 2 # Define known values. - self.si_sld = 2.073 - self.sio2_sld = 3.41 + self.si_sld = 2.073 + self.sio2_sld = 3.41 self.popc_hg_vol = 320.9 self.popc_tg_vol = 881.64 - self.popc_hg_sl = neutron_scattering_length('C10H18NO8P').real - self.popc_tg_sl = neutron_scattering_length('C32H64').real - self.water_vol = 30.4 + self.popc_hg_sl = neutron_scattering_length('C10H18NO8P').real + self.popc_tg_sl = neutron_scattering_length('C32H64').real + self.water_vol = 30.4 # Calculate the SLD of the tails. self.tg_sld = (self.popc_tg_sl / self.popc_tg_vol) * 1e6 # Define the varying parameters of the model. - self.si_rough = refnx.analysis.Parameter(2, 'Si/SiO2 Roughness', (1,8)) - self.sio2_thick = refnx.analysis.Parameter(14.7, 'SiO2 Thickness', (5,20)) - self.sio2_rough = refnx.analysis.Parameter(2, 'SiO2/POPC Roughness', (1,8)) - self.sio2_solv = refnx.analysis.Parameter(0.245, 'SiO2 Hydration', (0,1)) - self.popc_apm = refnx.analysis.Parameter(49.9, 'POPC Area Per Molecule', (30,60)) - self.bilayer_rough = refnx.analysis.Parameter(6.57, 'Bilayer Roughness', (1,8)) - self.bilayer_solv = refnx.analysis.Parameter(0.074, 'Bilayer Hydration', (0,1)) - self.hg_waters = refnx.analysis.Parameter(3.59, 'Headgroup Bound Waters', (0,20)) - - self.params = [self.si_rough, - self.sio2_thick, - self.sio2_rough, - self.sio2_solv, - self.popc_apm, - self.bilayer_rough, - self.bilayer_solv, - self.hg_waters] - - self.underlayer_params = [self.si_rough, - self.sio2_thick, - self.popc_apm, - self.bilayer_rough, - self.bilayer_solv, - self.hg_waters] + self.si_rough = refnx.analysis.Parameter( + 2, 'Si/SiO2 Roughness', (1, 8) + ) + self.sio2_thick = refnx.analysis.Parameter( + 14.7, 'SiO2 Thickness', (5, 20) + ) + self.sio2_rough = refnx.analysis.Parameter( + 2, 'SiO2/POPC Roughness', (1, 8) + ) + self.sio2_solv = refnx.analysis.Parameter( + 0.245, 'SiO2 Hydration', (0, 1) + ) + self.popc_apm = refnx.analysis.Parameter( + 49.9, 'POPC Area Per Molecule', (30, 60) + ) + self.bilayer_rough = refnx.analysis.Parameter( + 6.57, 'Bilayer Roughness', (1, 8) + ) + self.bilayer_solv = refnx.analysis.Parameter( + 0.074, 'Bilayer Hydration', (0, 1) + ) + self.hg_waters = refnx.analysis.Parameter( + 3.59, 'Headgroup Bound Waters', (0, 20) + ) + + self.params = [ + self.si_rough, + self.sio2_thick, + self.sio2_rough, + self.sio2_solv, + self.popc_apm, + self.bilayer_rough, + self.bilayer_solv, + self.hg_waters, + ] + + self.underlayer_params = [ + self.si_rough, + self.sio2_thick, + self.popc_apm, + self.bilayer_rough, + self.bilayer_solv, + self.hg_waters, + ] # Vary all of the parameters defined above. for param in self.params: - param.vary=True + param.vary = True # Call the BaseLipid constructor. super().__init__() @@ -130,9 +151,9 @@ def __init__(self): def _create_objectives(self): """Creates objectives corresponding to each measured contrast.""" # Define scattering lengths and densities of D2O and H2O. - d2o_sl = 2e-4 + d2o_sl = 2e-4 d2o_sld = 6.19 - h2o_sl = -1.64e-5 + h2o_sl = -1.64e-5 h2o_sld = -0.5227 D2O = refnx.reflect.SLD(d2o_sld) @@ -141,7 +162,7 @@ def _create_objectives(self): # Relate headgroup bound waters to scattering lengths and volumes. hg_water_d2o_sl = self.hg_waters * d2o_sl hg_water_h2o_sl = self.hg_waters * h2o_sl - hg_water_vol = self.hg_waters * self.water_vol + hg_water_vol = self.hg_waters * self.water_vol # Add to the headgroup volumes and scattering lengths in both contrast. vol_hg = self.popc_hg_vol + hg_water_vol @@ -150,32 +171,53 @@ def _create_objectives(self): popc_hg_sl_h2o = self.popc_hg_sl + hg_water_h2o_sl # Calculate the SLD of the headgroup in both contrast cases - sld_hg_d2o = (popc_hg_sl_d2o / vol_hg) * 1e6 # SLD = sum b / v + sld_hg_d2o = (popc_hg_sl_d2o / vol_hg) * 1e6 # SLD = sum b / v sld_hg_h2o = (popc_hg_sl_h2o / vol_hg) * 1e6 # Calculate thickness from headgroup volume over lipid APM. - hg_thick = vol_hg / self.popc_apm # Thickness = v / APM + hg_thick = vol_hg / self.popc_apm # Thickness = v / APM # Calculate the thickness of the tailgroup tg_thick = self.popc_tg_vol / self.popc_apm # Define the layers of the structure. substrate = refnx.reflect.SLD(self.si_sld) - sio2 = refnx.reflect.Slab(self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=self.sio2_solv) - - inner_hg_d2o = refnx.reflect.Slab(hg_thick, sld_hg_d2o, self.sio2_rough, vfsolv=self.bilayer_solv) - outer_hg_d2o = refnx.reflect.Slab(hg_thick, sld_hg_d2o, self.bilayer_rough, vfsolv=self.bilayer_solv) - inner_hg_h2o = refnx.reflect.Slab(hg_thick, sld_hg_d2o, self.sio2_rough, vfsolv=self.bilayer_solv) - outer_hg_h2o = refnx.reflect.Slab(hg_thick, sld_hg_h2o, self.bilayer_rough, vfsolv=self.bilayer_solv) - - tg = refnx.reflect.Slab(tg_thick, self.tg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv) + sio2 = refnx.reflect.Slab( + self.sio2_thick, + self.sio2_sld, + self.si_rough, + vfsolv=self.sio2_solv, + ) + + inner_hg_d2o = refnx.reflect.Slab( + hg_thick, sld_hg_d2o, self.sio2_rough, vfsolv=self.bilayer_solv + ) + outer_hg_d2o = refnx.reflect.Slab( + hg_thick, sld_hg_d2o, self.bilayer_rough, vfsolv=self.bilayer_solv + ) + inner_hg_h2o = refnx.reflect.Slab( + hg_thick, sld_hg_d2o, self.sio2_rough, vfsolv=self.bilayer_solv + ) + outer_hg_h2o = refnx.reflect.Slab( + hg_thick, sld_hg_h2o, self.bilayer_rough, vfsolv=self.bilayer_solv + ) + + tg = refnx.reflect.Slab( + tg_thick, self.tg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv + ) # Structure corresponding to measuring the Si/D2O interface. si_D2O_structure = substrate | sio2 | D2O(rough=self.sio2_rough) # Two structures corresponding to each measured contrast. - si_POPC_D2O_structure = substrate | sio2 | inner_hg_d2o | tg | tg | outer_hg_d2o | D2O(rough=self.bilayer_rough) - si_POPC_H2O_structure = substrate | sio2 | inner_hg_h2o | tg | tg | outer_hg_h2o | H2O(rough=self.bilayer_rough) + si_POPC_D2O_structure = ( + substrate | sio2 | inner_hg_d2o | tg | tg + | outer_hg_d2o | D2O(rough=self.bilayer_rough) + ) + si_POPC_H2O_structure = ( + substrate | sio2 | inner_hg_h2o | tg | tg + | outer_hg_h2o | H2O(rough=self.bilayer_rough) + ) self.structures = [si_D2O_structure, si_POPC_D2O_structure, @@ -185,10 +227,9 @@ def _create_objectives(self): self.objectives = [] for i, structure in enumerate(self.structures): # Define the model. - model = refnx.reflect.ReflectModel(structure, - scale=self.scales[i], - bkg=self.bkgs[i], - dq=self.dq) + model = refnx.reflect.ReflectModel( + structure, scale=self.scales[i], bkg=self.bkgs[i], dq=self.dq + ) # Load the measured data. filename = '{}.dat'.format(self.labels[i]) file_path = os.path.join(self.data_path, filename) @@ -210,10 +251,10 @@ def _using_conditions(self, contrast_sld, underlayers=None): """ # Calculate the SLD of the headgroup with the given contrast SLD. - hg_sld = contrast_sld*0.27 + 1.98*0.73 + hg_sld = contrast_sld * 0.27 + 1.98 * 0.73 # Calculate the headgroup and tailgroup thicknesses. - vol_hg = self.popc_hg_vol + self.hg_waters*self.water_vol + vol_hg = self.popc_hg_vol + self.hg_waters * self.water_vol hg_thick = vol_hg / self.popc_apm tg_thick = self.popc_tg_vol / self.popc_apm @@ -221,22 +262,37 @@ def _using_conditions(self, contrast_sld, underlayers=None): substrate = refnx.reflect.SLD(self.si_sld) solution = refnx.reflect.SLD(contrast_sld)(rough=self.bilayer_rough) - inner_hg = refnx.reflect.Slab(hg_thick, hg_sld, self.sio2_rough, vfsolv=self.bilayer_solv) - outer_hg = refnx.reflect.Slab(hg_thick, hg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv) - tg = refnx.reflect.Slab(tg_thick, self.tg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv) + inner_hg = refnx.reflect.Slab( + hg_thick, hg_sld, self.sio2_rough, vfsolv=self.bilayer_solv + ) + outer_hg = refnx.reflect.Slab( + hg_thick, hg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv + ) + tg = refnx.reflect.Slab( + tg_thick, self.tg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv + ) # Add underlayers if specified. if underlayers is None: - sio2 = refnx.reflect.Slab(self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=self.sio2_solv) + sio2 = refnx.reflect.Slab( + self.sio2_thick, + self.sio2_sld, + self.si_rough, + vfsolv=self.sio2_solv, + ) return substrate | sio2 | inner_hg | tg | tg | outer_hg | solution else: # Use 0% hydration for the SiO2 layer. - sio2 = refnx.reflect.Slab(self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=0) + sio2 = refnx.reflect.Slab( + self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=0 + ) structure = substrate | sio2 # Add each underlayer with given thickness and SLD. for thick, sld in underlayers: - underlayer = refnx.reflect.SLD(sld)(thick, 2) # Default 2 roughness. + underlayer = refnx.reflect.SLD(sld)( + thick, 2 + ) # Default 2 roughness. structure |= underlayer return structure | inner_hg | tg | tg | outer_hg | solution @@ -276,12 +332,12 @@ class BilayerDMPC(BaseLipid): objectives (list): objectives corresponding to each measured contrast. """ + def __init__(self): self.name = 'DMPC_bilayer' - self.data_path = os.path.join(os.path.dirname(__file__), - '..', - 'data', - 'DMPC_bilayer') + self.data_path = os.path.join( + os.path.dirname(__file__), '..', 'data', 'DMPC_bilayer' + ) self.labels = ['Si-D2O', 'Si-DMPC-D2O', 'Si-DMPC-H2O'] self.distances = np.linspace(-20, 95, 500) @@ -290,46 +346,66 @@ def __init__(self): self.dq = 2 # Define known values. - self.si_sld = 2.073 - self.sio2_sld = 3.41 + self.si_sld = 2.073 + self.sio2_sld = 3.41 self.dmpc_hg_vol = 320.9 self.dmpc_tg_vol = 783.3 - self.dmpc_hg_sl = 6.41e-4 - self.dmpc_tg_sl = -3.08e-4 - self.water_vol = 30.4 + self.dmpc_hg_sl = 6.41e-4 + self.dmpc_tg_sl = -3.08e-4 + self.water_vol = 30.4 # Calculate the SLD of the tails. self.tg_sld = (self.dmpc_tg_sl / self.dmpc_tg_vol) * 1e6 # Define the varying parameters of the model. - self.si_rough = refnx.analysis.Parameter(2, 'Si/SiO2 Roughness', (1,8)) - self.sio2_thick = refnx.analysis.Parameter(14.7, 'SiO2 Thickness', (5,20)) - self.sio2_rough = refnx.analysis.Parameter(2, 'SiO2/DMPC Roughness', (1,8)) - self.sio2_solv = refnx.analysis.Parameter(0.245, 'SiO2 Hydration', (0,1)) - self.dmpc_apm = refnx.analysis.Parameter(49.9, 'DMPC Area Per Molecule', (30,60)) - self.bilayer_rough = refnx.analysis.Parameter(6.57, 'Bilayer Roughness', (1,8)) - self.bilayer_solv = refnx.analysis.Parameter(0.074, 'Bilayer Hydration', (0,1)) - self.hg_waters = refnx.analysis.Parameter(3.59, 'Headgroup Bound Waters', (0,20)) - - self.params = [self.si_rough, - self.sio2_thick, - self.sio2_rough, - self.sio2_solv, - self.dmpc_apm, - self.bilayer_rough, - self.bilayer_solv, - self.hg_waters] - - self.underlayer_params = [self.si_rough, - self.sio2_thick, - self.dmpc_apm, - self.bilayer_rough, - self.bilayer_solv, - self.hg_waters] + self.si_rough = refnx.analysis.Parameter( + 2, 'Si/SiO2 Roughness', (1, 8) + ) + self.sio2_thick = refnx.analysis.Parameter( + 14.7, 'SiO2 Thickness', (5, 20) + ) + self.sio2_rough = refnx.analysis.Parameter( + 2, 'SiO2/DMPC Roughness', (1, 8) + ) + self.sio2_solv = refnx.analysis.Parameter( + 0.245, 'SiO2 Hydration', (0, 1) + ) + self.dmpc_apm = refnx.analysis.Parameter( + 49.9, 'DMPC Area Per Molecule', (30, 60) + ) + self.bilayer_rough = refnx.analysis.Parameter( + 6.57, 'Bilayer Roughness', (1, 8) + ) + self.bilayer_solv = refnx.analysis.Parameter( + 0.074, 'Bilayer Hydration', (0, 1) + ) + self.hg_waters = refnx.analysis.Parameter( + 3.59, 'Headgroup Bound Waters', (0, 20) + ) + + self.params = [ + self.si_rough, + self.sio2_thick, + self.sio2_rough, + self.sio2_solv, + self.dmpc_apm, + self.bilayer_rough, + self.bilayer_solv, + self.hg_waters, + ] + + self.underlayer_params = [ + self.si_rough, + self.sio2_thick, + self.dmpc_apm, + self.bilayer_rough, + self.bilayer_solv, + self.hg_waters, + ] # Vary all of the parameters defined above. for param in self.params: - param.vary=True + param.vary = True # Call the BaseLipid constructor. super().__init__() @@ -337,9 +413,9 @@ def __init__(self): def _create_objectives(self): """Creates objectives corresponding to each measured contrast.""" # Define scattering lengths and densities of D2O and H2O. - d2o_sl = 2e-4 + d2o_sl = 2e-4 d2o_sld = 6.19 - h2o_sl = -1.64e-5 + h2o_sl = -1.64e-5 h2o_sld = -0.5227 D2O = refnx.reflect.SLD(d2o_sld) @@ -348,7 +424,7 @@ def _create_objectives(self): # Relate headgroup bound waters to scattering lengths and volumes. hg_water_d2o_sl = self.hg_waters * d2o_sl hg_water_h2o_sl = self.hg_waters * h2o_sl - hg_water_vol = self.hg_waters * self.water_vol + hg_water_vol = self.hg_waters * self.water_vol # Add to the headgroup volumes and scattering lengths in both contrast. vol_hg = self.dmpc_hg_vol + hg_water_vol @@ -357,45 +433,77 @@ def _create_objectives(self): dmpc_hg_sl_h2o = self.dmpc_hg_sl + hg_water_h2o_sl # Calculate the SLD of the headgroup in both contrast cases - sld_hg_d2o = (dmpc_hg_sl_d2o / vol_hg) * 1e6 # SLD = sum b / v + sld_hg_d2o = (dmpc_hg_sl_d2o / vol_hg) * 1e6 # SLD = sum b / v sld_hg_h2o = (dmpc_hg_sl_h2o / vol_hg) * 1e6 # Calculate thickness from headgroup volume over lipid APM. - hg_thick = vol_hg / self.dmpc_apm # Thickness = v / APM + hg_thick = vol_hg / self.dmpc_apm # Thickness = v / APM # Calculate the thickness of the tailgroup tg_thick = self.dmpc_tg_vol / self.dmpc_apm # Define the layers of the structure. substrate = refnx.reflect.SLD(self.si_sld) - sio2 = refnx.reflect.Slab(self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=self.sio2_solv) - - inner_hg_d2o = refnx.reflect.Slab(hg_thick, sld_hg_d2o, self.sio2_rough, vfsolv=self.bilayer_solv) - outer_hg_d2o = refnx.reflect.Slab(hg_thick, sld_hg_d2o, self.bilayer_rough, vfsolv=self.bilayer_solv) - inner_hg_h2o = refnx.reflect.Slab(hg_thick, sld_hg_d2o, self.sio2_rough, vfsolv=self.bilayer_solv) - outer_hg_h2o = refnx.reflect.Slab(hg_thick, sld_hg_h2o, self.bilayer_rough, vfsolv=self.bilayer_solv) - - tg = refnx.reflect.Slab(tg_thick, self.tg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv) + sio2 = refnx.reflect.Slab( + self.sio2_thick, + self.sio2_sld, + self.si_rough, + vfsolv=self.sio2_solv, + ) + + inner_hg_d2o = refnx.reflect.Slab( + hg_thick, sld_hg_d2o, self.sio2_rough, vfsolv=self.bilayer_solv + ) + outer_hg_d2o = refnx.reflect.Slab( + hg_thick, sld_hg_d2o, self.bilayer_rough, vfsolv=self.bilayer_solv + ) + inner_hg_h2o = refnx.reflect.Slab( + hg_thick, sld_hg_d2o, self.sio2_rough, vfsolv=self.bilayer_solv + ) + outer_hg_h2o = refnx.reflect.Slab( + hg_thick, sld_hg_h2o, self.bilayer_rough, vfsolv=self.bilayer_solv + ) + + tg = refnx.reflect.Slab( + tg_thick, self.tg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv + ) # Structure corresponding to measuring the Si/D2O interface. si_D2O_structure = substrate | sio2 | D2O(rough=self.sio2_rough) # Two structures corresponding to each measured contrast. - si_DMPC_D2O_structure = substrate | sio2 | inner_hg_d2o | tg | tg | outer_hg_d2o | D2O(rough=self.bilayer_rough) - si_DMPC_H2O_structure = substrate | sio2 | inner_hg_h2o | tg | tg | outer_hg_h2o | H2O(rough=self.bilayer_rough) - - self.structures = [si_D2O_structure, - si_DMPC_D2O_structure, - si_DMPC_H2O_structure] + si_DMPC_D2O_structure = ( + substrate + | sio2 + | inner_hg_d2o + | tg + | tg + | outer_hg_d2o + | D2O(rough=self.bilayer_rough) + ) + si_DMPC_H2O_structure = ( + substrate + | sio2 + | inner_hg_h2o + | tg + | tg + | outer_hg_h2o + | H2O(rough=self.bilayer_rough) + ) + + self.structures = [ + si_D2O_structure, + si_DMPC_D2O_structure, + si_DMPC_H2O_structure, + ] # Iterate over each measured structure. self.objectives = [] for i, structure in enumerate(self.structures): # Define the model. - model = refnx.reflect.ReflectModel(structure, - scale=self.scales[i], - bkg=self.bkgs[i], - dq=self.dq) + model = refnx.reflect.ReflectModel( + structure, scale=self.scales[i], bkg=self.bkgs[i], dq=self.dq + ) # Load the measured data. filename = '{}.dat'.format(self.labels[i]) file_path = os.path.join(self.data_path, filename) @@ -417,10 +525,10 @@ def _using_conditions(self, contrast_sld, underlayers=None): """ # Calculate the SLD of the headgroup with the given contrast SLD. - hg_sld = contrast_sld*0.27 + 1.98*0.73 + hg_sld = contrast_sld * 0.27 + 1.98 * 0.73 # Calculate the headgroup and tailgroup thicknesses. - vol_hg = self.dmpc_hg_vol + self.hg_waters*self.water_vol + vol_hg = self.dmpc_hg_vol + self.hg_waters * self.water_vol hg_thick = vol_hg / self.dmpc_apm tg_thick = self.dmpc_tg_vol / self.dmpc_apm @@ -428,26 +536,42 @@ def _using_conditions(self, contrast_sld, underlayers=None): substrate = refnx.reflect.SLD(self.si_sld) solution = refnx.reflect.SLD(contrast_sld)(rough=self.bilayer_rough) - inner_hg = refnx.reflect.Slab(hg_thick, hg_sld, self.sio2_rough, vfsolv=self.bilayer_solv) - outer_hg = refnx.reflect.Slab(hg_thick, hg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv) - tg = refnx.reflect.Slab(tg_thick, self.tg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv) + inner_hg = refnx.reflect.Slab( + hg_thick, hg_sld, self.sio2_rough, vfsolv=self.bilayer_solv + ) + outer_hg = refnx.reflect.Slab( + hg_thick, hg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv + ) + tg = refnx.reflect.Slab( + tg_thick, self.tg_sld, self.bilayer_rough, vfsolv=self.bilayer_solv + ) # Add underlayers if specified. if underlayers is None: - sio2 = refnx.reflect.Slab(self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=self.sio2_solv) + sio2 = refnx.reflect.Slab( + self.sio2_thick, + self.sio2_sld, + self.si_rough, + vfsolv=self.sio2_solv, + ) return substrate | sio2 | inner_hg | tg | tg | outer_hg | solution else: # Use 0% hydration for the SiO2 layer. - sio2 = refnx.reflect.Slab(self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=0) + sio2 = refnx.reflect.Slab( + self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=0 + ) structure = substrate | sio2 # Add each underlayer with given thickness and SLD. for thick, sld in underlayers: - underlayer = refnx.reflect.SLD(sld)(thick, 2) # Default 2 roughness. + underlayer = refnx.reflect.SLD(sld)( + thick, 2 + ) # Default 2 roughness. structure |= underlayer return structure | inner_hg | tg | tg | outer_hg | solution + class BilayerDPPC(BaseLipid): """Defines a model describing a 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) and Ra lipopolysaccharide (LPS) bilayer defined by a single @@ -460,7 +584,8 @@ class BilayerDPPC(BaseLipid): distances (numpy.ndarray): SLD profile x-axis range. contrast_slds (list): SLD of each measured contrast. scale (float): experimental scale factor for measured contrasts. - bkgs (list): level of instrument background noise for each measured contrast. + bkgs (list): level of instrument background noise for each measured + contrast. dq (float): instrument resolution. si_sld (float): SLD of silicon substrate. sio2_sld (float): SLD of silicon oxide. @@ -488,12 +613,12 @@ class BilayerDPPC(BaseLipid): objectives (list): objectives corresponding to each measured contrast. """ + def __init__(self): self.name = 'DPPC_RaLPS_bilayer' - self.data_path = os.path.join(os.path.dirname(__file__), - '..', - 'data', - 'DPPC_RaLPS_bilayer') + self.data_path = os.path.join( + os.path.dirname(__file__), '..', 'data', 'DPPC_RaLPS_bilayer' + ) self.labels = ['dDPPC-RaLPS-D2O', 'dDPPC-RaLPS-SMW', 'dDPPC-RaLPS-H2O'] self.distances = np.linspace(-30, 110, 500) @@ -513,48 +638,78 @@ def __init__(self): self.core_h2o_sld = 2.01 # Define the varying parameters of the model. - self.si_rough = refnx.analysis.Parameter(5.5, 'Si/SiO2 Roughness', (3,8)) - self.sio2_thick = refnx.analysis.Parameter(13.4, 'SiO2 Thickness', (10,30)) - self.sio2_rough = refnx.analysis.Parameter(3.2, 'SiO2/Bilayer Roughness', (2,5)) - self.sio2_solv = refnx.analysis.Parameter(0.038, 'SiO2 Hydration', (0,0.5)) - self.inner_hg_thick = refnx.analysis.Parameter(9.0, 'Inner Headgroup Thickness', (5,20)) - self.inner_hg_solv = refnx.analysis.Parameter(0.39, 'Inner Headgroup Hydration', (0,1)) - self.bilayer_rough = refnx.analysis.Parameter(4.0, 'Bilayer Roughness', (0,12)) - self.inner_tg_thick = refnx.analysis.Parameter(16.7, 'Inner Tailgroup Thickness', (10,20)) - self.outer_tg_thick = refnx.analysis.Parameter(14.9, 'Outer Tailgroup Thickness', (10,20)) - self.tg_solv = refnx.analysis.Parameter(0.0085, 'Tailgroup Hydration', (0,1)) - self.core_thick = refnx.analysis.Parameter(28.7, 'Core Thickness', (0,50)) - self.core_solv = refnx.analysis.Parameter(0.26, 'Core Hydration', (0,1)) - self.asym_value = refnx.analysis.Parameter(0.95, 'Asymmetry Value', (0,1)) - - self.params = [self.si_rough, - self.sio2_thick, - self.sio2_rough, - self.sio2_solv, - self.inner_hg_thick, - self.inner_hg_solv, - self.bilayer_rough, - self.inner_tg_thick, - self.outer_tg_thick, - self.tg_solv, - self.core_thick, - self.core_solv, - self.asym_value] - - self.underlayer_params = [self.si_rough, - self.sio2_thick, - self.inner_hg_thick, - self.inner_hg_solv, - self.bilayer_rough, - self.inner_tg_thick, - self.outer_tg_thick, - self.tg_solv, - self.core_thick, - self.core_solv] + self.si_rough = refnx.analysis.Parameter( + 5.5, 'Si/SiO2 Roughness', (3, 8) + ) + self.sio2_thick = refnx.analysis.Parameter( + 13.4, 'SiO2 Thickness', (10, 30) + ) + self.sio2_rough = refnx.analysis.Parameter( + 3.2, 'SiO2/Bilayer Roughness', (2, 5) + ) + self.sio2_solv = refnx.analysis.Parameter( + 0.038, 'SiO2 Hydration', (0, 0.5) + ) + self.inner_hg_thick = refnx.analysis.Parameter( + 9.0, 'Inner Headgroup Thickness', (5, 20) + ) + self.inner_hg_solv = refnx.analysis.Parameter( + 0.39, 'Inner Headgroup Hydration', (0, 1) + ) + self.bilayer_rough = refnx.analysis.Parameter( + 4.0, 'Bilayer Roughness', (0, 12) + ) + self.inner_tg_thick = refnx.analysis.Parameter( + 16.7, 'Inner Tailgroup Thickness', (10, 20) + ) + self.outer_tg_thick = refnx.analysis.Parameter( + 14.9, 'Outer Tailgroup Thickness', (10, 20) + ) + self.tg_solv = refnx.analysis.Parameter( + 0.0085, 'Tailgroup Hydration', (0, 1) + ) + self.core_thick = refnx.analysis.Parameter( + 28.7, 'Core Thickness', (0, 50) + ) + self.core_solv = refnx.analysis.Parameter( + 0.26, 'Core Hydration', (0, 1) + ) + self.asym_value = refnx.analysis.Parameter( + 0.95, 'Asymmetry Value', (0, 1) + ) + + self.params = [ + self.si_rough, + self.sio2_thick, + self.sio2_rough, + self.sio2_solv, + self.inner_hg_thick, + self.inner_hg_solv, + self.bilayer_rough, + self.inner_tg_thick, + self.outer_tg_thick, + self.tg_solv, + self.core_thick, + self.core_solv, + self.asym_value, + ] + + self.underlayer_params = [ + self.si_rough, + self.sio2_thick, + self.inner_hg_thick, + self.inner_hg_solv, + self.bilayer_rough, + self.inner_tg_thick, + self.outer_tg_thick, + self.tg_solv, + self.core_thick, + self.core_solv, + ] # Vary all of the parameters defined above. for param in self.params: - param.vary=True + param.vary = True # Call the BaseLipid constructor. super().__init__() @@ -562,17 +717,18 @@ def __init__(self): def _create_objectives(self): """Creates objectives corresponding to each measured contrast.""" # Define structures for each contrast. - self.structures = [self._using_conditions(contrast_sld) - for contrast_sld in self.contrast_slds] + self.structures = [ + self._using_conditions(contrast_sld) + for contrast_sld in self.contrast_slds + ] # Iterate over each measured structure. self.objectives = [] for i, structure in enumerate(self.structures): # Define the model. - model = refnx.reflect.ReflectModel(structure, - scale=self.scale, - bkg=self.bkgs[i], - dq=self.dq) + model = refnx.reflect.ReflectModel( + structure, scale=self.scale, bkg=self.bkgs[i], dq=self.dq + ) # Load the measured data. filename = '{}.dat'.format(self.labels[i]) file_path = os.path.join(self.data_path, filename) @@ -595,46 +751,93 @@ def _using_conditions(self, contrast_sld, underlayers=None): """ # Calculate core SLD with the given contrast SLD. contrast_point = (contrast_sld + 0.56) / (6.35 + 0.56) - core_sld = contrast_point*self.core_d2o_sld + (1-contrast_point)*self.core_h2o_sld + core_sld = ( + contrast_point * self.core_d2o_sld + + (1 - contrast_point) * self.core_h2o_sld + ) # Use asymmetry to define inner and outer tailgroup SLDs. - inner_tg_sld = self.asym_value*self.dppc_tg_sld + (1-self.asym_value)*self.lps_tg_sld - outer_tg_sld = self.asym_value*self.lps_tg_sld + (1-self.asym_value)*self.dppc_tg_sld + inner_tg_sld = ( + self.asym_value * self.dppc_tg_sld + + (1 - self.asym_value) * self.lps_tg_sld + ) + outer_tg_sld = ( + self.asym_value * self.lps_tg_sld + + (1 - self.asym_value) * self.dppc_tg_sld + ) # Define the layers of the model. substrate = refnx.reflect.SLD(self.si_sld) - inner_hg = refnx.reflect.Slab(self.inner_hg_thick, self.dppc_hg_sld, self.sio2_rough, vfsolv=self.inner_hg_solv) - inner_tg = refnx.reflect.Slab(self.inner_tg_thick, inner_tg_sld, self.bilayer_rough, vfsolv=self.tg_solv) - outer_tg = refnx.reflect.Slab(self.outer_tg_thick, outer_tg_sld, self.bilayer_rough, vfsolv=self.tg_solv) - - core = refnx.reflect.Slab(self.core_thick, core_sld, self.bilayer_rough, vfsolv=self.core_solv) + inner_hg = refnx.reflect.Slab( + self.inner_hg_thick, + self.dppc_hg_sld, + self.sio2_rough, + vfsolv=self.inner_hg_solv, + ) + inner_tg = refnx.reflect.Slab( + self.inner_tg_thick, + inner_tg_sld, + self.bilayer_rough, + vfsolv=self.tg_solv, + ) + outer_tg = refnx.reflect.Slab( + self.outer_tg_thick, + outer_tg_sld, + self.bilayer_rough, + vfsolv=self.tg_solv, + ) + + core = refnx.reflect.Slab( + self.core_thick, + core_sld, + self.bilayer_rough, + vfsolv=self.core_solv, + ) solution = refnx.reflect.SLD(contrast_sld)(0, self.bilayer_rough) # Add the underlayers if specified. if underlayers is None: - sio2 = refnx.reflect.Slab(self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=self.sio2_solv) - return substrate | sio2 | inner_hg | inner_tg | outer_tg | core | solution + sio2 = refnx.reflect.Slab( + self.sio2_thick, + self.sio2_sld, + self.si_rough, + vfsolv=self.sio2_solv, + ) + return ( + substrate + | sio2 + | inner_hg + | inner_tg + | outer_tg + | core + | solution + ) else: # Use 0% hydration for the SiO2 layer. - sio2 = refnx.reflect.Slab(self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=0) + sio2 = refnx.reflect.Slab( + self.sio2_thick, self.sio2_sld, self.si_rough, vfsolv=0 + ) structure = substrate | sio2 # Add each underlayer with given thickness and SLD. for thick, sld in underlayers: - underlayer = refnx.reflect.SLD(sld)(thick, 2) # Default 2 roughness. + underlayer = refnx.reflect.SLD(sld)( + thick, 2 + ) # Default 2 roughness. structure |= underlayer return structure | inner_hg | inner_tg | outer_tg | core | solution + if __name__ == '__main__': save_path = './results' # Save the SLD and reflectivity profiles of the POPC bilayer. - #dmpc_bilayer = BilayerPOPC() - #dmpc_bilayer.sld_profile(save_path) - #dmpc_bilayer.reflectivity_profile(save_path) + # dmpc_bilayer = BilayerPOPC() + # dmpc_bilayer.sld_profile(save_path) + # dmpc_bilayer.reflectivity_profile(save_path) # Save the SLD and reflectivity profiles of the DMPC bilayer. dmpc_bilayer = BilayerDMPC() diff --git a/hogben/models/magnetic.py b/hogben/models/magnetic.py index 1f84324a..cb2e97db 100644 --- a/hogben/models/magnetic.py +++ b/hogben/models/magnetic.py @@ -1,5 +1,4 @@ import os -import sys import matplotlib.pyplot as plt @@ -39,8 +38,10 @@ class SampleYIG(BaseSample, VariableUnderlayer): pt_rough (bumps.parameter.Parameter): air/platinum roughness. pt_mag (bumps.parameter.Parameter): platinum layer magnetic SLD. intermediary_sld (bumps.parameter.Parameter): intermediary layer SLD. - intermediary_thick (bumps.parameter.Parameter): intermediary layer thickness. - intermediary_rough (bumps.parameter.Parameter): platinum/intermediary roughness. + intermediary_thick (bumps.parameter.Parameter): intermediary layer + thickness. + intermediary_rough (bumps.parameter.Parameter): platinum/intermediary + roughness. yig_sld (bumps.parameter.Parameter): YIG layer SLD. yig_thick (bumps.parameter.Parameter): YIG layer thickness. yig_rough (bumps.parameter.Parameter): intermediary/YIG roughness. @@ -49,15 +50,16 @@ class SampleYIG(BaseSample, VariableUnderlayer): yag_rough (bumps.parameter.Parameter): YIG/YAG roughness. params (list): parameters of the model. structure (refl1d.model.Stack): Refl1D representation of sample. - experiment (refl1d.experiment.Experiment): fittable experiment for sample. + experiment (refl1d.experiment.Experiment): fittable experiment for + sample. """ + def __init__(self): self.name = 'YIG_sample' - self.data_path = os.path.join(os.path.dirname(__file__), - '..', - 'data', - 'YIG_sample') + self.data_path = os.path.join( + os.path.dirname(__file__), '..', 'data', 'YIG_sample' + ) self.labels = ['Up', 'Down'] self.mag_angle = 90 @@ -68,20 +70,34 @@ def __init__(self): # Define the parameters of the model. self.pt_sld = bumps.parameter.Parameter(5.646, name='Pt SLD') self.pt_thick = bumps.parameter.Parameter(21.08, name='Pt Thickness') - self.pt_rough = bumps.parameter.Parameter(8.211, name='Air/Pt Roughness') + self.pt_rough = bumps.parameter.Parameter( + 8.211, name='Air/Pt Roughness' + ) self.pt_mag = bumps.parameter.Parameter(0, name='Pt Magnetic SLD') - self.intermediary_sld = bumps.parameter.Parameter(4.678, name='Intermediary SLD') - self.intermediary_thick = bumps.parameter.Parameter(19.67, name='Intermediary Thickness') - self.intermediary_rough = bumps.parameter.Parameter(2, name='Pt/Intermediary Roughness') + self.intermediary_sld = bumps.parameter.Parameter( + 4.678, name='Intermediary SLD' + ) + self.intermediary_thick = bumps.parameter.Parameter( + 19.67, name='Intermediary Thickness' + ) + self.intermediary_rough = bumps.parameter.Parameter( + 2, name='Pt/Intermediary Roughness' + ) self.yig_sld = bumps.parameter.Parameter(5.836, name='YIG SLD') self.yig_thick = bumps.parameter.Parameter(713.8, name='YIG Thickness') - self.yig_rough = bumps.parameter.Parameter(13.55, name='Intermediary/YIG Roughness') - self.yig_mag = bumps.parameter.Parameter(0.349, name='YIG Magnetic SLD') + self.yig_rough = bumps.parameter.Parameter( + 13.55, name='Intermediary/YIG Roughness' + ) + self.yig_mag = bumps.parameter.Parameter( + 0.349, name='YIG Magnetic SLD' + ) self.yag_sld = bumps.parameter.Parameter(5.304, name='YAG SLD') - self.yag_rough = bumps.parameter.Parameter(30, name='YIG/YAG Roughness') + self.yag_rough = bumps.parameter.Parameter( + 30, name='YIG/YAG Roughness' + ) """ self.params = [self.pt_sld, @@ -127,11 +143,13 @@ def __create_experiment(self): file_path_up = os.path.join(self.data_path, 'YIG_up.dat') file_path_down = os.path.join(self.data_path, 'YIG_down.dat') - pp = refl1d.probe.load4(file_path_up, sep='\t', - intensity=self.scale, background=self.bkg) + pp = refl1d.probe.load4( + file_path_up, sep='\t', intensity=self.scale, background=self.bkg + ) - mm = refl1d.probe.load4(file_path_down, sep='\t', - intensity=self.scale, background=self.bkg) + mm = refl1d.probe.load4( + file_path_down, sep='\t', intensity=self.scale, background=self.bkg + ) # Set the resolution to be constant dQ/q. self.__set_dq(pp) @@ -142,7 +160,9 @@ def __create_experiment(self): # Define the experiment using the probe and sample structure. self.structure = self.using_conditions() - self.experiment = refl1d.experiment.Experiment(sample=self.structure, probe=probe) + self.experiment = refl1d.experiment.Experiment( + sample=self.structure, probe=probe + ) def using_conditions(self, yig_thick=None, pt_thick=None): """Creates a structure representing the YIG sample measured using @@ -169,18 +189,32 @@ def using_conditions(self, yig_thick=None, pt_thick=None): assert pt_thick >= 0 # 1 uB/atom = 1.638 - pt_magnetism = refl1d.magnetism.Magnetism(rhoM=self.pt_mag, thetaM=self.mag_angle) - yig_magnetism = refl1d.magnetism.Magnetism(rhoM=self.yig_mag, thetaM=self.mag_angle) + pt_magnetism = refl1d.magnetism.Magnetism( + rhoM=self.pt_mag, thetaM=self.mag_angle + ) + yig_magnetism = refl1d.magnetism.Magnetism( + rhoM=self.yig_mag, thetaM=self.mag_angle + ) # Define the sample structure. air = refl1d.material.SLD(rho=0, name='Air') - pt = refl1d.material.SLD(rho=self.pt_sld, name='Pt')(self.pt_thick, self.pt_rough, magnetism=pt_magnetism) - intermediary = refl1d.material.SLD(rho=self.intermediary_sld, name='Intermediary')(self.intermediary_thick, self.intermediary_rough) - yig = refl1d.material.SLD(rho=self.yig_sld, name='YIG')(yig_thick, self.yig_rough, magnetism=yig_magnetism) - yag = refl1d.material.SLD(rho=self.yag_sld, name='Substrate')(0, self.yag_rough) + pt = refl1d.material.SLD(rho=self.pt_sld, name='Pt')( + self.pt_thick, self.pt_rough, magnetism=pt_magnetism + ) + intermediary = refl1d.material.SLD( + rho=self.intermediary_sld, name='Intermediary' + )(self.intermediary_thick, self.intermediary_rough) + yig = refl1d.material.SLD(rho=self.yig_sld, name='YIG')( + yig_thick, self.yig_rough, magnetism=yig_magnetism + ) + yag = refl1d.material.SLD(rho=self.yag_sld, name='Substrate')( + 0, self.yag_rough + ) # Define the added platinum thickness, if applicable. - pt_extra = refl1d.material.SLD(rho=self.pt_sld, name='Pt Extra')(pt_thick, 0) + pt_extra = refl1d.material.SLD(rho=self.pt_sld, name='Pt Extra')( + pt_thick, 0 + ) # Add the extra platinum layer if requested. if pt_thick == 0: @@ -207,8 +241,8 @@ def angle_info(self, angle_times, contrasts=None): mp=False, mm=True) # Calculate the Fisher information matrix. - qs = [data[:,0] for data in datasets] - counts = [data[:,3] for data in datasets] + qs = [data[:, 0] for data in datasets] + counts = [data[:, 3] for data in datasets] return fisher(qs, self.params, counts, models) def underlayer_info(self, angle_times, yig_thick, pt_thick): @@ -234,8 +268,8 @@ def underlayer_info(self, angle_times, yig_thick, pt_thick): mp=False, mm=True) # Calculate the Fisher information matrix. - qs = [data[:,0] for data in datasets] - counts = [data[:,3] for data in datasets] + qs = [data[:, 0] for data in datasets] + counts = [data[:, 3] for data in datasets] return fisher(qs, self.params, counts, models) def __set_dq(self, probe): @@ -246,7 +280,7 @@ def __set_dq(self, probe): """ # Transform the resolution from refnx to Refl1D format. - dq = self.dq / (100*np.sqrt(8*np.log(2))) + dq = self.dq / (100 * np.sqrt(8 * np.log(2))) q_array = probe.Q @@ -256,9 +290,9 @@ def __set_dq(self, probe): # Adjust probe calculation for constant resolution. argmin, argmax = np.argmin(q_array), np.argmax(q_array) - probe.calc_Qo = np.linspace(q_array[argmin] - 3.5*dq_array[argmin], - q_array[argmax] + 3.5*dq_array[argmax], - 21*len(q_array)) + probe.calc_Qo = np.linspace(q_array[argmin] - 3.5 * dq_array[argmin], + q_array[argmax] + 3.5 * dq_array[argmax], + 21 * len(q_array)) def sld_profile(self, save_path): """Plots the SLD profile of the sample. @@ -270,7 +304,7 @@ def sld_profile(self, save_path): # Get the SLD profile values. z, slds, _, slds_mag, _ = self.experiment.magnetic_smooth_profile() - fig = plt.figure(figsize=(8,6)) + fig = plt.figure(figsize=(8, 6)) ax = fig.add_subplot(111) # Plot the SLD profile. @@ -295,7 +329,7 @@ def reflectivity_profile(self, save_path): save_path (str): path to directory to save profile to. """ - fig = plt.figure(figsize=(6,7)) + fig = plt.figure(figsize=(6, 7)) ax = fig.add_subplot(111) colours = ['b', 'g'] @@ -310,12 +344,12 @@ def reflectivity_profile(self, save_path): ax.errorbar(probe.Q, probe.R, probe.dR, marker='o', ms=2, lw=0, elinewidth=0.5, capsize=0.5, color=colours[count], - label=self.labels[count]+' Data') + label=self.labels[count] + ' Data') # Plot the model reflectivity. ax.plot(probe.Q, qr[1], zorder=20, color=colours[count], - label=self.labels[count]+' Fitted') + label=self.labels[count] + ' Fitted') count += 1 @@ -350,7 +384,8 @@ def nested_sampling(self, angle_times, save_path, filename, dynamic=False): # Save the sampling corner plot. save_path = os.path.join(save_path, self.name) - save_plot(fig, save_path, 'nested_sampling_'+filename) + save_plot(fig, save_path, 'nested_sampling_' + filename) + if __name__ == '__main__': save_path = '../results' diff --git a/hogben/models/monolayers.py b/hogben/models/monolayers.py index 8efc3ab1..62501ed3 100644 --- a/hogben/models/monolayers.py +++ b/hogben/models/monolayers.py @@ -1,4 +1,5 @@ import os +from typing import Optional import matplotlib.pyplot as plt @@ -37,15 +38,15 @@ class MonolayerDPPG(BaseLipid): params (list): varying model parameters. """ - def __init__(self, deuterated=False): + + def __init__(self, deuterated: bool = False): self.name = 'DPPG_monolayer' - self.data_path = os.path.join(os.path.dirname(__file__), - '..', - 'data', - 'DPPG_monolayer') + self.data_path = os.path.join( + os.path.dirname(__file__), '..', 'data', 'DPPG_monolayer' + ) self.labels = ['hDPPG-D2O', 'dDPPG-NRW', 'hDPPG-NRW'] - self.distances = None #np.linspace(-25, 90, 500) + self.distances = None # np.linspace(-25, 90, 500) self.deuterated = deuterated self.scales = [1.8899, 1.8832, 1.8574] @@ -53,47 +54,68 @@ def __init__(self, deuterated=False): self.dq = 3 # Define the varying parameters of the model. - self.air_tg_rough = refnx.analysis.Parameter( 5.0000, 'Air/Tailgroup Roughness', ( 5, 8)) - self.lipid_apm = refnx.analysis.Parameter(54.1039, 'Lipid Area Per Molecule', (30, 80)) - self.hg_waters = refnx.analysis.Parameter( 6.6874, 'Headgroup Bound Waters', ( 0, 20)) - self.monolayer_rough = refnx.analysis.Parameter( 2.0233, 'Monolayer Roughness', ( 0, 10)) - self.lipid_vfsolv = refnx.analysis.Parameter( 0.9046, 'Lipid Hydration', ( 0, 1)) - self.protein_thick = refnx.analysis.Parameter(32.6858, 'Protein Thickness', ( 0, 80)) - self.protein_vfsolv = refnx.analysis.Parameter( 0.5501, 'Protein Hydration', ( 0, 100)) - self.water_rough = refnx.analysis.Parameter( 3.4590, 'Protein/Water Roughness', ( 0, 15)) + self.air_tg_rough = refnx.analysis.Parameter( + 5.0000, 'Air/Tailgroup Roughness', (5, 8) + ) + self.lipid_apm = refnx.analysis.Parameter( + 54.1039, 'Lipid Area Per Molecule', (30, 80) + ) + self.hg_waters = refnx.analysis.Parameter( + 6.6874, 'Headgroup Bound Waters', (0, 20) + ) + self.monolayer_rough = refnx.analysis.Parameter( + 2.0233, 'Monolayer Roughness', (0, 10) + ) + self.lipid_vfsolv = refnx.analysis.Parameter( + 0.9046, 'Lipid Hydration', (0, 1) + ) + self.protein_thick = refnx.analysis.Parameter( + 32.6858, 'Protein Thickness', (0, 80) + ) + self.protein_vfsolv = refnx.analysis.Parameter( + 0.5501, 'Protein Hydration', (0, 100) + ) + self.water_rough = refnx.analysis.Parameter( + 3.4590, 'Protein/Water Roughness', (0, 15) + ) # We are only interested in the lipid APM. self.params = [self.lipid_apm] # Vary all of the parameters defined above. for param in self.params: - param.vary=True + param.vary = True # Call the BaseLipid constructor. super().__init__() - def _create_objectives(self, protein=True): + def _create_objectives(self, protein: bool = True) -> None: """Creates objectives corresponding to each measured contrast. Args: - protein (bool): whether the protein is included or not in the model. + protein (bool): whether the protein is included or not in the + model. """ # SLDs of null-reflecting water (NRW) and D2O. nrw, d2o = 0.1, 6.35 # Define the measured structures. - self.structures = [self._using_conditions(d2o, deuterated=False, protein=protein), - self._using_conditions(nrw, deuterated=True, protein=protein), - self._using_conditions(nrw, deuterated=False, protein=protein)] + self.structures = [ + self._using_conditions(d2o, deuterated=False, protein=protein), + self._using_conditions(nrw, deuterated=True, protein=protein), + self._using_conditions(nrw, deuterated=False, protein=protein), + ] self.objectives = [] for i, structure in enumerate(self.structures): # Define the model. - model = refnx.reflect.ReflectModel(structure, - scale=self.scales[i], - bkg=self.bkgs[i]*self.scales[i], - dq=self.dq) + model = refnx.reflect.ReflectModel( + structure, + scale=self.scales[i], + bkg=self.bkgs[i] * self.scales[i], + dq=self.dq, + ) # Load the measured data. filename = '{}.dat'.format(self.labels[i]) file_path = os.path.join(self.data_path, filename) @@ -102,8 +124,12 @@ def _create_objectives(self, protein=True): # Combine model and data into an objective that can be fitted. self.objectives.append(refnx.analysis.Objective(model, data)) - def _using_conditions(self, contrast_sld, underlayers=None, - deuterated=None, protein=False): + def _using_conditions(self, + contrast_sld: float, + underlayers: Optional[list] = None, + deuterated: Optional[bool] = None, + protein: bool = False + ) -> refnx.reflect.Structure: """Creates a structure representing the monolayer measured using given measurement conditions. @@ -128,7 +154,7 @@ def _using_conditions(self, contrast_sld, underlayers=None, contrast_sld *= 1e-6 # Define known SLDs of D2O and H2O - d2o_sld = 6.35e-6 + d2o_sld = 6.35e-6 h2o_sld = -0.56e-6 # Define known protein SLDs for D2O and H2O. @@ -136,49 +162,49 @@ def _using_conditions(self, contrast_sld, underlayers=None, protein_h2o_sld = 1.9e-6 # Define neutron scattering lengths. - carbon_sl = 0.6646e-4 - oxygen_sl = 0.5843e-4 - hydrogen_sl = -0.3739e-4 - phosphorus_sl = 0.5130e-4 - deuterium_sl = 0.6671e-4 + carbon_sl = 0.6646e-4 + oxygen_sl = 0.5843e-4 + hydrogen_sl = -0.3739e-4 + phosphorus_sl = 0.5130e-4 + deuterium_sl = 0.6671e-4 # Calculate the total scattering length in each fragment. - ch2_sl = 1*carbon_sl + 2*hydrogen_sl - ch3_sl = 1*carbon_sl + 3*hydrogen_sl - co2_sl = 1*carbon_sl + 2*oxygen_sl - c3h5_sl = 3*carbon_sl + 5*hydrogen_sl - po4_sl = 1*phosphorus_sl + 4*oxygen_sl - h2o_sl = 2*hydrogen_sl + 1*oxygen_sl - d2o_sl = 2*deuterium_sl + 1*oxygen_sl - cd2_sl = 1*carbon_sl + 2*deuterium_sl - cd3_sl = 1*carbon_sl + 3*deuterium_sl + ch2_sl = 1 * carbon_sl + 2 * hydrogen_sl + ch3_sl = 1 * carbon_sl + 3 * hydrogen_sl + co2_sl = 1 * carbon_sl + 2 * oxygen_sl + c3h5_sl = 3 * carbon_sl + 5 * hydrogen_sl + po4_sl = 1 * phosphorus_sl + 4 * oxygen_sl + h2o_sl = 2 * hydrogen_sl + 1 * oxygen_sl + d2o_sl = 2 * deuterium_sl + 1 * oxygen_sl + cd2_sl = 1 * carbon_sl + 2 * deuterium_sl + cd3_sl = 1 * carbon_sl + 3 * deuterium_sl # Volumes of each fragment. - ch2_vol = 28.1 - ch3_vol = 52.7/2 - co2_vol = 39.0 - c3h5_vol = 68.8 - po4_vol = 53.7 + ch2_vol = 28.1 + ch3_vol = 52.7 / 2 + co2_vol = 39.0 + c3h5_vol = 68.8 + po4_vol = 53.7 water_vol = 30.4 # Calculate volumes from components. - hg_vol = po4_vol + 2*c3h5_vol + 2*co2_vol - tg_vol = 28*ch2_vol + 2*ch3_vol + hg_vol = po4_vol + 2 * c3h5_vol + 2 * co2_vol + tg_vol = 28 * ch2_vol + 2 * ch3_vol # Calculate mole fraction of D2O from the bulk SLD. x = (contrast_sld - h2o_sld) / (d2o_sld - h2o_sld) # Calculate 'average' scattering length sum per water molecule. - sl_sum_water = x*d2o_sl + (1-x)*h2o_sl + sl_sum_water = x * d2o_sl + (1 - x) * h2o_sl # Calculate scattering length sums for the other fragments. - sl_sum_hg = po4_sl + 2*c3h5_sl + 2*co2_sl - sl_sum_tg_h = 28*ch2_sl + 2*ch3_sl - sl_sum_tg_d = 28*cd2_sl + 2*cd3_sl + sl_sum_hg = po4_sl + 2 * c3h5_sl + 2 * co2_sl + sl_sum_tg_h = 28 * ch2_sl + 2 * ch3_sl + sl_sum_tg_d = 28 * cd2_sl + 2 * cd3_sl # Need to include number of hydrating water molecules in headgroup. - hg_vol += water_vol*self.hg_waters - sl_sum_hg += sl_sum_water*self.hg_waters + hg_vol += water_vol * self.hg_waters + sl_sum_hg += sl_sum_water * self.hg_waters hg_sld = sl_sum_hg / hg_vol @@ -189,13 +215,18 @@ def _using_conditions(self, contrast_sld, underlayers=None, # If including the protein in the model. if protein: # Calculate the SLD of the protein. - protein_sld = x*protein_d2o_sld + (1-x)*protein_h2o_sld - protein = refnx.reflect.SLD(protein_sld*1e6, name='Protein')(self.protein_thick, self.monolayer_rough, self.protein_vfsolv) + protein_sld = x * protein_d2o_sld + (1 - x) * protein_h2o_sld + protein = refnx.reflect.SLD(protein_sld * 1e6, name='Protein')( + self.protein_thick, self.monolayer_rough, self.protein_vfsolv + ) # Adjust headgroup and tailgroup SLDs by lipid hydration. - hg_sld = self.lipid_vfsolv*hg_sld + (1-self.lipid_vfsolv)*protein_sld - tg_h_sld = self.lipid_vfsolv*tg_h_sld + (1-self.lipid_vfsolv)*protein_sld - tg_d_sld = self.lipid_vfsolv*tg_d_sld + (1-self.lipid_vfsolv)*protein_sld + hg_sld = (self.lipid_vfsolv * hg_sld + + (1 - self.lipid_vfsolv) * protein_sld) + tg_h_sld = (self.lipid_vfsolv * tg_h_sld + + (1 - self.lipid_vfsolv) * protein_sld) + tg_d_sld = (self.lipid_vfsolv * tg_d_sld + + (1 - self.lipid_vfsolv) * protein_sld) # Tailgroup and headgroup thicknesses. tg_thick = tg_vol / self.lipid_apm @@ -203,10 +234,18 @@ def _using_conditions(self, contrast_sld, underlayers=None, # Define the structure. air = refnx.reflect.SLD(0, name='Air') - tg_h = refnx.reflect.SLD(tg_h_sld*1e6, name='Hydrogenated Tailgroup')(tg_thick, self.air_tg_rough) - tg_d = refnx.reflect.SLD(tg_d_sld*1e6, name='Deuterated Tailgroup')(tg_thick, self.air_tg_rough) - hg = refnx.reflect.SLD(hg_sld*1e6, name='Headgroup')(hg_thick, self.monolayer_rough) - water = refnx.reflect.SLD(contrast_sld*1e6, name='Water')(0, self.water_rough) + tg_h = refnx.reflect.SLD( + tg_h_sld * 1e6, name='Hydrogenated Tailgroup' + )(tg_thick, self.air_tg_rough) + tg_d = refnx.reflect.SLD(tg_d_sld * 1e6, name='Deuterated Tailgroup')( + tg_thick, self.air_tg_rough + ) + hg = refnx.reflect.SLD(hg_sld * 1e6, name='Headgroup')( + hg_thick, self.monolayer_rough + ) + water = refnx.reflect.SLD(contrast_sld * 1e6, name='Water')( + 0, self.water_rough + ) # Add either hydrogenated or deuterated tailgroups. structure = air @@ -221,7 +260,7 @@ def _using_conditions(self, contrast_sld, underlayers=None, else: return structure | hg | water - def sld_profile(self, save_path): + def sld_profile(self, save_path: str) -> None: """Plots the SLD profiles of the monolayer sample. Args: @@ -230,12 +269,14 @@ def sld_profile(self, save_path): """ # Plot the SLD profile without the protein. self._create_objectives(protein=False) - super().sld_profile(save_path, 'sld_profile_no_protein', ylim=(-0.6, 7.5)) + super().sld_profile(save_path, 'sld_profile_no_protein', + ylim=(-0.6, 7.5)) # Plot the SLD profile with the protein. self._create_objectives(protein=True) super().sld_profile(save_path, 'sld_profile_protein', ylim=(-0.6, 7.5)) + if __name__ == '__main__': save_path = '../results' diff --git a/hogben/models/parsing.py b/hogben/models/parsing.py index cd3898cb..b40eef25 100644 --- a/hogben/models/parsing.py +++ b/hogben/models/parsing.py @@ -38,7 +38,8 @@ def _fuse(mol1: dict, mol2: dict, w: int = 1) -> dict: def _parse(formula: str) -> Tuple[dict, int]: """ :param formula: Chemical formula as a string - :return: Tuple containing; formula as a dictwith occurences of each atom and an iterator. + :return: Tuple containing; formula as a dictwith occurences of each atom + and an iterator. """ q = [] mol = {} @@ -75,7 +76,7 @@ def _parse(formula: str) -> Tuple[dict, int]: def parse_formula(formula: str) -> dict: - """ + """ :param formula: Chemical formula as a string :return: Formula as a dict with occurences of each atom. """ diff --git a/hogben/models/samples.py b/hogben/models/samples.py index 1e8d13e7..d427f58e 100644 --- a/hogben/models/samples.py +++ b/hogben/models/samples.py @@ -1,6 +1,4 @@ import os -import sys - import matplotlib.pyplot as plt import numpy as np @@ -22,7 +20,7 @@ from hogben.utils import fisher, Sampler, save_plot from hogben.models.base import BaseSample -plt.rcParams['figure.figsize'] = (9,7) +plt.rcParams['figure.figsize'] = (9, 7) plt.rcParams['figure.dpi'] = 600 @@ -36,6 +34,7 @@ class Sample(BaseSample): params (list): varying parameters of sample. """ + def __init__(self, structure): self.structure = structure self.name = structure.name @@ -60,12 +59,18 @@ def __vary_structure(structure, bound_size=0.2): # Vary the SLD and thickness of each component (layer). for component in structure[1:-1]: sld = component.sld.real - sld_bounds = (sld.value*(1-bound_size), sld.value*(1+bound_size)) + sld_bounds = ( + sld.value * (1 - bound_size), + sld.value * (1 + bound_size), + ) sld.setp(vary=True, bounds=sld_bounds) params.append(sld) thick = component.thick - thick_bounds = (thick.value*(1-bound_size), thick.value*(1+bound_size)) + thick_bounds = ( + thick.value * (1 - bound_size), + thick.value * (1 + bound_size), + ) thick.setp(vary=True, bounds=thick_bounds) params.append(thick) @@ -74,11 +79,11 @@ def __vary_structure(structure, bound_size=0.2): # Vary the SLD and thickness of each component (layer). for component in structure[1:-1]: sld = component.material.rho - sld.pmp(bound_size*100) + sld.pmp(bound_size * 100) params.append(sld) thick = component.thickness - thick.pmp(bound_size*100) + thick.pmp(bound_size * 100) params.append(thick) # Otherwise, the structure is invalid. @@ -100,7 +105,7 @@ def angle_info(self, angle_times, contrasts=None): """ # Return the Fisher information matrix calculated from simulated data. model, data = simulate(self.structure, angle_times) - qs, counts, models = [data[:,0]], [data[:,3]], [model] + qs, counts, models = [data[:, 0]], [data[:, 3]], [model] return fisher(qs, self.params, counts, models) def sld_profile(self, save_path): @@ -116,8 +121,8 @@ def sld_profile(self, save_path): # Determine if the structure was defined in Refl1D. elif isinstance(self.structure, refl1d.model.Stack): - q = np.geomspace(0.005, 0.3, 500) # This is not used. - scale, bkg, dq = 1, 1e-6, 2 # These are not used. + q = np.geomspace(0.005, 0.3, 500) # This is not used. + scale, bkg, dq = 1, 1e-6, 2 # These are not used. experiment = refl1d_experiment(self.structure, q, scale, bkg, dq) z, slds, _ = experiment.smooth_profile() @@ -141,8 +146,16 @@ def sld_profile(self, save_path): save_path = os.path.join(save_path, self.name) save_plot(fig, save_path, 'sld_profile') - def reflectivity_profile(self, save_path, q_min=0.005, q_max=0.4, - points=500, scale=1, bkg=1e-7, dq=2): + def reflectivity_profile( + self, + save_path: str, + q_min: float = 0.005, + q_max: float = 0.4, + points: int = 500, + scale: float = 1, + bkg: float = 1e-7, + dq: float = 2, + ) -> None: """Plots the reflectivity profile of the sample. Args: @@ -160,7 +173,9 @@ def reflectivity_profile(self, save_path, q_min=0.005, q_max=0.4, # Determine if the structure was defined in refnx. if isinstance(self.structure, refnx.reflect.Structure): - model = refnx.reflect.ReflectModel(self.structure, scale=scale, bkg=bkg, dq=dq) + model = refnx.reflect.ReflectModel( + self.structure, scale=scale, bkg=bkg, dq=dq + ) # Determine if the structure was defined in Refl1D. elif isinstance(self.structure, refl1d.model.Stack): @@ -189,7 +204,11 @@ def reflectivity_profile(self, save_path, q_min=0.005, q_max=0.4, save_path = os.path.join(save_path, self.name) save_plot(fig, save_path, 'reflectivity_profile') - def nested_sampling(self, angle_times, save_path, filename, dynamic=False): + def nested_sampling(self, + angle_times: list, + save_path: str, + filename: str, + dynamic: bool = False) -> None: """Runs nested sampling on simulated data of the sample. Args: @@ -204,7 +223,9 @@ def nested_sampling(self, angle_times, save_path, filename, dynamic=False): # Determine if the structure was defined in refnx. if isinstance(self.structure, refnx.reflect.Structure): - dataset = refnx.reflect.ReflectDataset([data[:,0], data[:,1], data[:,2]]) + dataset = refnx.reflect.ReflectDataset( + [data[:, 0], data[:, 1], data[:, 2]] + ) objective = refnx.anaylsis.Objective(model, dataset) # Determine if the structure was defined in Refl1D. @@ -221,7 +242,7 @@ def nested_sampling(self, angle_times, save_path, filename, dynamic=False): # Save the sampling corner plot. save_path = os.path.join(save_path, self.name) - save_plot(fig, save_path, filename+'_nested_sampling') + save_plot(fig, save_path, filename + '_nested_sampling') def to_refl1d(self): """Converts the refnx structure to an equivalent Refl1D structure.""" @@ -234,7 +255,8 @@ def to_refl1d(self): name, sld = component.name, component.sld.real.value, thick, rough = component.thick.value, component.rough.value - # Add the component in the opposite direction to the refnx definition. + # Add the component in the opposite direction to the refnx + # definition. layer = refl1d.material.SLD(rho=sld, name=name)(thick, rough) structure = layer | structure @@ -253,13 +275,15 @@ def to_refnx(self): name, sld = component.name, component.material.rho.value, thick, rough = component.thickness.value, component.interface.value - # Add the component in the opposite direction to the Refl1D definition. + # Add the component in the opposite direction to the Refl1D + # definition. structure |= refnx.reflect.SLD(sld, name=name)(thick, rough) # Update the current structure to use the new version. structure.name = self.structure.name self.structure = structure + def simple_sample(): """Defines a 2-layer simple sample. @@ -276,6 +300,7 @@ def simple_sample(): structure.name = 'simple_sample' return Sample(structure) + def many_param_sample(): """Defines a 5-layer sample with many parameters. @@ -295,6 +320,7 @@ def many_param_sample(): structure.name = 'many_param_sample' return Sample(structure) + def thin_layer_sample_1(): """Defines a 2-layer sample with thin layers. @@ -311,6 +337,7 @@ def thin_layer_sample_1(): structure.name = 'thin_layer_sample_1' return Sample(structure) + def thin_layer_sample_2(): """Defines a 3-layer sample with thin layers. @@ -328,6 +355,7 @@ def thin_layer_sample_2(): structure.name = 'thin_layer_sample_2' return Sample(structure) + def similar_sld_sample_1(): """Defines a 2-layer sample with layers of similar SLD. @@ -344,6 +372,7 @@ def similar_sld_sample_1(): structure.name = 'similar_sld_sample_1' return Sample(structure) + def similar_sld_sample_2(): """Defines a 3-layer sample with layers of similar SLD. @@ -361,6 +390,7 @@ def similar_sld_sample_2(): structure.name = 'similar_sld_sample_2' return Sample(structure) + if __name__ == '__main__': save_path = '../results' diff --git a/hogben/optimise.py b/hogben/optimise.py index 7b48ca7b..719fb12d 100644 --- a/hogben/optimise.py +++ b/hogben/optimise.py @@ -1,11 +1,12 @@ -import os -import sys - import numpy as np from scipy.optimize import differential_evolution, NonlinearConstraint -from hogben.models.base import VariableAngle, VariableContrast, VariableUnderlayer +from hogben.models.base import ( + VariableAngle, + VariableContrast, + VariableUnderlayer, +) class Optimiser: @@ -15,12 +16,20 @@ class Optimiser: sample (base.BaseSample): sample to optimise an experiment for. """ + def __init__(self, sample): self.sample = sample - def optimise_angle_times(self, num_angles, contrasts=[], total_time=1000, - angle_bounds=(0.2, 4), points=100, - workers=-1, verbose=True): + def optimise_angle_times( + self, + num_angles: int, + contrasts: list = [], + total_time: float = 1000, + angle_bounds: tuple = (0.2, 4), + points: int = 100, + workers: int = -1, + verbose: bool = True, + ) -> tuple: """Optimises the measurement angles and associated counting times of an experiment, given a fixed time budget. @@ -42,15 +51,19 @@ def optimise_angle_times(self, num_angles, contrasts=[], total_time=1000, assert isinstance(self.sample, VariableAngle) # Define bounds on each condition to optimise (angles and time splits). - bounds = [angle_bounds]*num_angles + [(0, 1)]*num_angles + bounds = [angle_bounds] * num_angles + [(0, 1)] * num_angles # Arguments for the optimisation function. args = [num_angles, contrasts, points, total_time] # Constrain the counting times to sum to the fixed time budget. # Also constrain the angles to be in non-decreasing order. - sum_of_splits = lambda x: sum(x[num_angles:]) - non_decreasing = lambda x: int(np.all(np.diff(x[:num_angles]) >= 0)) + def sum_of_splits(x): + return sum(x[num_angles:]) + + def non_decreasing(x): + return int(np.all(np.diff(x[:num_angles]) >= 0)) + constraints = [NonlinearConstraint(sum_of_splits, 1, 1), NonlinearConstraint(non_decreasing, 1, 1)] @@ -59,9 +72,15 @@ def optimise_angle_times(self, num_angles, contrasts=[], total_time=1000, constraints, args, workers, verbose) return res[:num_angles], res[num_angles:], val - def optimise_contrasts(self, num_contrasts, angle_splits, - total_time=1000, contrast_bounds=(-0.56, 6.36), - workers=-1, verbose=True): + def optimise_contrasts( + self, + num_contrasts: int, + angle_splits: list, + total_time: float = 1000, + contrast_bounds: tuple = (-0.56, 6.36), + workers: int = -1, + verbose: bool = True, + ) -> tuple: """Finds the optimal contrasts, given a fixed time budget. Args: @@ -82,26 +101,40 @@ def optimise_contrasts(self, num_contrasts, angle_splits, # Define the bounds on each condition to optimise # (contrast SLDs and time splits). - bounds = [contrast_bounds]*num_contrasts + [(0, 1)]*num_contrasts + bounds = [contrast_bounds] * num_contrasts + [(0, 1)] * num_contrasts # Constrain the counting times to sum to the fixed time budget. # Also constrain the contrasts to be in non-decreasing order. - sum_of_splits = lambda x: sum(x[num_contrasts:]) - non_decreasing = lambda x: int(np.all(np.diff(x[:num_contrasts]) >= 0)) - constraints = [NonlinearConstraint(sum_of_splits, 1, 1), - NonlinearConstraint(non_decreasing, 1, 1)] + def sum_of_splits(x): + return sum(x[num_contrasts:]) + + def non_decreasing(x): + return int(np.all(np.diff(x[:num_contrasts]) >= 0)) + + constraints = [ + NonlinearConstraint(sum_of_splits, 1, 1), + NonlinearConstraint(non_decreasing, 1, 1), + ] # Arguments for the optimisation function. args = [num_contrasts, angle_splits, total_time] # Optimise contrasts and counting time splits, and return the results. - res, val = Optimiser.__optimise(self._contrasts_func, bounds, - constraints, args, workers, verbose) + res, val = Optimiser.__optimise( + self._contrasts_func, bounds, constraints, args, workers, verbose + ) return res[:num_contrasts], res[num_contrasts:], val - def optimise_underlayers(self, num_underlayers, angle_times, contrasts, - thick_bounds=(0, 500), sld_bounds=(1, 9), - workers=-1, verbose=True): + def optimise_underlayers( + self, + num_underlayers, + angle_times, + contrasts, + thick_bounds=(0, 500), + sld_bounds=(1, 9), + workers=-1, + verbose=True, + ) -> tuple: """Finds the optimal underlayer thicknesses and SLDs of a sample. Args: @@ -123,17 +156,25 @@ def optimise_underlayers(self, num_underlayers, angle_times, contrasts, # Define bounds on each condition to optimise # (underlayer thicknesses and SLDs). - bounds = [thick_bounds]*num_underlayers + [sld_bounds]*num_underlayers + bounds = [thick_bounds] * num_underlayers + [ + sld_bounds + ] * num_underlayers # Arguments for the optimisation function. args = [num_underlayers, angle_times, contrasts] # Optimise underlayer thicknesses and SLDs, and return the results. - res, val = Optimiser.__optimise(self._underlayers_func, bounds, [], - args, workers, verbose) + res, val = Optimiser.__optimise( + self._underlayers_func, bounds, [], args, workers, verbose + ) return res[:num_underlayers], res[num_underlayers:], val - def _angle_times_func(self, x, num_angles, contrasts, points, total_time): + def _angle_times_func(self, + x: list, + num_angles: int, + contrasts: list, + points: int, + total_time: float) -> float: """Defines the function for optimising an experiment's measurement angles and associated counting times. @@ -149,8 +190,10 @@ def _angle_times_func(self, x, num_angles, contrasts, points, total_time): """ # Extract the angles and counting times from given list, `x`. - angle_times = [(x[i], points, total_time*x[num_angles+i]) - for i in range(num_angles)] + angle_times = [ + (x[i], points, total_time * x[num_angles + i]) + for i in range(num_angles) + ] # Calculate the Fisher information matrix. g = self.sample.angle_info(angle_times, contrasts) @@ -158,7 +201,11 @@ def _angle_times_func(self, x, num_angles, contrasts, points, total_time): # Return negative of the minimum eigenvalue as algorithm is minimising. return -np.linalg.eigvalsh(g)[0] - def _contrasts_func(self, x, num_contrasts, angle_splits, total_time): + def _contrasts_func(self, + x: list, + num_contrasts: int, + angle_splits: type, + total_time: float) -> float: """Defines the function for optimising an experiment's contrasts. Args: @@ -178,8 +225,10 @@ def _contrasts_func(self, x, num_contrasts, angle_splits, total_time): # Iterate over each contrast. for i in range(num_contrasts): # Calculate proportion of the total counting time for each angle. - angle_times = [(angle, points, total_time*x[num_contrasts+i]*split) - for angle, points, split in angle_splits] + angle_times = [ + (angle, points, total_time * x[num_contrasts + i] * split) + for angle, points, split in angle_splits + ] # Add to the initial Fisher information matrix. g += self.sample.contrast_info(angle_times, [x[i]]) @@ -187,7 +236,11 @@ def _contrasts_func(self, x, num_contrasts, angle_splits, total_time): # Return negative of the minimum eigenvalue as algorithm is minimising. return -np.linalg.eigvalsh(g)[0] - def _underlayers_func(self, x, num_underlayers, angle_times, contrasts): + def _underlayers_func(self, + x: list, + num_underlayers: int, + angle_times: type, + contrasts: list) -> float: """Defines the function for optimising an experiment's underlayers. Args: @@ -201,8 +254,9 @@ def _underlayers_func(self, x, num_underlayers, angle_times, contrasts): """ # Extract the underlayer thicknesses and SLDs from the given `x` list. - underlayers = [(x[i], x[num_underlayers+i]) - for i in range(num_underlayers)] + underlayers = [ + (x[i], x[num_underlayers + i]) for i in range(num_underlayers) + ] # Calculate the Fisher information matrix using the conditions. g = self.sample.underlayer_info(angle_times, contrasts, underlayers) @@ -211,7 +265,12 @@ def _underlayers_func(self, x, num_underlayers, angle_times, contrasts): return -np.linalg.eigvalsh(g)[0] @staticmethod - def __optimise(func, bounds, constraints, args, workers, verbose): + def __optimise(func: callable, + bounds: list, + constraints: list, + args: list, + workers: int, + verbose: bool) -> tuple: """Optimises a given `func` using the differential evolution global optimisation algorithm. @@ -229,7 +288,8 @@ def __optimise(func, bounds, constraints, args, workers, verbose): """ # Run differential evolution on the given optimisation function. res = differential_evolution(func, bounds, constraints=constraints, - args=args, polish=False, tol=0.001, - updating='deferred', workers=workers, - disp=verbose) + args=args, polish=False, tol=0.001, + updating='deferred', workers=workers, + disp=verbose) + return res.x, res.fun diff --git a/hogben/simulate.py b/hogben/simulate.py index b9c7db63..5ac5e2d8 100644 --- a/hogben/simulate.py +++ b/hogben/simulate.py @@ -27,15 +27,17 @@ def direct_beam_path(inst_or_path: str = 'OFFSPEC', beam file or the local path """ - non_pol_instr = {'OFFSPEC': 'OFFSPEC_non_polarised_old.dat', - 'SURF': 'SURF_non_polarised.dat', - 'POLREF': 'POLREF_non_polarised.dat', - 'INTER': 'INTER_non_polarised.dat' - } - - pol_instr = {'OFFSPEC': 'OFFSPEC_polarised_old.dat', - 'POLREF': 'POLREF_polarised.dat' - } + non_pol_instr = { + 'OFFSPEC': 'OFFSPEC_non_polarised_old.dat', + 'SURF': 'SURF_non_polarised.dat', + 'POLREF': 'POLREF_non_polarised.dat', + 'INTER': 'INTER_non_polarised.dat', + } + + pol_instr = { + 'OFFSPEC': 'OFFSPEC_polarised_old.dat', + 'POLREF': 'POLREF_polarised.dat', + } # Check if the key isn't in the dictionary and assume # a local filepath @@ -43,26 +45,34 @@ def direct_beam_path(inst_or_path: str = 'OFFSPEC', if os.path.isfile(inst_or_path): return inst_or_path else: - msg = "Please provide an instrument name or a valid local filepath" + msg = 'Please provide an instrument name or a valid local filepath' raise FileNotFoundError(str(msg)) path = importlib_resources.files('hogben.data.directbeams').joinpath( - non_pol_instr[inst_or_path]) + non_pol_instr[inst_or_path] + ) if polarised: - path = importlib_resources.files('hogben.data.directbeams' - ).joinpath(pol_instr[inst_or_path]) + path = importlib_resources.files('hogben.data.directbeams').joinpath( + pol_instr[inst_or_path] + ) return path -def simulate_magnetic(sample: Union['refnx.reflect.Stucture', - 'refl1d.model.Stack'], - angle_times: np.ndarray, scale: float = 1.0, - bkg: float = 5e-7, dq: float = 2, mm: bool = True, - mp: bool = True, pm: bool = True, pp: bool = True, - angle_scale: float = 0.3, - inst_or_path: str = 'OFFSPEC') -> tuple: +def simulate_magnetic( + sample: Union['refnx.reflect.Stucture', 'refl1d.model.Stack'], + angle_times: np.ndarray, + scale: float = 1.0, + bkg: float = 5e-7, + dq: float = 2, + mm: bool = True, + mp: bool = True, + pm: bool = True, + pp: bool = True, + angle_scale: float = 0.3, + inst_or_path: str = 'OFFSPEC', +) -> tuple: """Simulates an experiment of a given magnetic `sample` measured over a number of angles. @@ -178,8 +188,8 @@ def simulate(sample: Union['refnx.reflect.Stucture', 'refl1d.model.Stack'], data[:, 2] = np.concatenate(dr) data[:, 3] = np.concatenate(counts) - data = data[(data != 0).all(1)] # Remove points of zero reflectivity. - data = data[data[:,0].argsort()] # Sort by Q. + data = data[(data != 0).all(1)] # Remove points of zero reflectivity. + data = data[data[:, 0].argsort()] # Sort by Q. # If there is no data after removing zeros, return no model. if len(data) == 0: @@ -212,11 +222,18 @@ def simulate(sample: Union['refnx.reflect.Stucture', 'refl1d.model.Stack'], return model, data -def _run_experiment(sample: Union['refnx.reflect.Stucture', - 'refl1d.model.Stack'], - angle: float, points: int, time: float, - scale: float, bkg: float, dq: float, directbeam_path: str, - angle_scale: float, spin_state: int) -> tuple: +def _run_experiment( + sample: Union['refnx.reflect.Stucture', 'refl1d.model.Stack'], + angle: float, + points: int, + time: float, + scale: float, + bkg: float, + dq: float, + directbeam_path: str, + angle_scale: float, + spin_state: int, +) -> tuple: """Simulates a single angle measurement of a given `sample`. Args: @@ -241,18 +258,19 @@ def _run_experiment(sample: Union['refnx.reflect.Stucture', wavelengths = direct_beam[:, 0] # 1st column is wavelength, 2nd is flux. # Adjust flux by measurement angle. - direct_flux = direct_beam[:, 1] * pow(angle/angle_scale, 2) + direct_flux = direct_beam[:, 1] * pow(angle / angle_scale, 2) # Calculate Q values. - q = 4*np.pi*np.sin(np.radians(angle)) / wavelengths + q = 4 * np.pi * np.sin(np.radians(angle)) / wavelengths # Bin Q's' in equally geometrically-spaced bins using flux as weighting. - q_bin_edges = np.geomspace(min(q), max(q), points+1) + q_bin_edges = np.geomspace(min(q), max(q), points + 1) flux_binned, _ = np.histogram(q, q_bin_edges, weights=direct_flux) # Get the bin centres. - q_binned = np.asarray([(q_bin_edges[i] + q_bin_edges[i+1]) / 2 - for i in range(points)]) + q_binned = np.asarray( + [(q_bin_edges[i] + q_bin_edges[i + 1]) / 2 for i in range(points)] + ) # Calculate the model reflectivity at each Q point. if isinstance(sample, refnx.reflect.Structure): @@ -275,7 +293,8 @@ def _run_experiment(sample: Union['refnx.reflect.Stucture', # Get the measured reflected count for each bin. # r_model accounts for background. - counts_reflected = np.random.poisson(r_model*counts_incident).astype(float) + counts_reflected = np.random.poisson(r_model * counts_incident).astype( + float) # Convert from count space to reflectivity space. # Point has zero reflectivity if there is no flux. @@ -334,10 +353,13 @@ def reflectivity(q: np.ndarray, model: Union['refnx.reflect.ReflectModel', return experiment.reflectivity()[1] -def refl1d_experiment(sample: 'refl1d.model.stack', q_array: np.ndarray, - scale: float, bkg: float, dq: float, - spin_state: Optional[int] = None)\ - -> refl1d.experiment.Experiment: +def refl1d_experiment(sample: 'refl1d.model.stack', + q_array: np.ndarray, + scale: float, + bkg: float, + dq: float, + spin_state: Optional[int] = None, + ) -> refl1d.experiment.Experiment: """Creates a Refl1D experiment for a given `sample` and `q_array`. Args: @@ -353,7 +375,7 @@ def refl1d_experiment(sample: 'refl1d.model.stack', q_array: np.ndarray, """ # Transform the resolution from refnx to Refl1D format. - refl1d_dq = dq / (100*np.sqrt(8*np.log(2))) + refl1d_dq = dq / (100 * np.sqrt(8 * np.log(2))) # Calculate the dq array and use it to define a QProbe. dq_array = q_array * refl1d_dq @@ -363,13 +385,15 @@ def refl1d_experiment(sample: 'refl1d.model.stack', q_array: np.ndarray, # Adjust probe calculation for constant dQ/Q resolution. argmin, argmax = np.argmin(q_array), np.argmax(q_array) - probe.calc_Qo = np.linspace(q_array[argmin] - 3.5*dq_array[argmin], - q_array[argmax] + 3.5*dq_array[argmax], - 21*len(q_array)) + probe.calc_Qo = np.linspace( + q_array[argmin] - 3.5 * dq_array[argmin], + q_array[argmax] + 3.5 * dq_array[argmax], + 21 * len(q_array), + ) # If the sample is magnetic, create a polarised QProbe. if sample.ismagnetic: - probes = [None]*4 + probes = [None] * 4 probes[spin_state] = probe probe = refl1d.probe.PolarizedQProbe(xs=probes, name='') probe.spin_state = spin_state diff --git a/hogben/tests/test_fisher.py b/hogben/tests/test_fisher.py index a7f14ab4..339664a8 100644 --- a/hogben/tests/test_fisher.py +++ b/hogben/tests/test_fisher.py @@ -66,13 +66,20 @@ def mock_refnx_model(): bounds """ # Parameters described as tuples: (value, lower bound, upper bound) - parameter_values = [(20, 15, 25), (50, 45, 55), (10, 7.5, 8.5), - (2, 1.5, 2.5)] + parameter_values = [ + (20, 15, 25), + (50, 45, 55), + (10, 7.5, 8.5), + (2, 1.5, 2.5), + ] # Fill parameter values and bounds parameters = [ - Mock(spec=refnx.analysis.Parameter, value=value, - bounds=Mock(lb=lb, ub=ub)) + Mock( + spec=refnx.analysis.Parameter, + value=value, + bounds=Mock(lb=lb, ub=ub), + ) for value, lb, ub in parameter_values ] model = Mock(spec=refnx.reflect.ReflectModel, xi=parameters) @@ -87,13 +94,20 @@ def mock_refl1d_model(): bounds """ # Parameters described as tuples: (value, lower bound, upper bound) - parameter_values = [(20, 15, 25), (50, 45, 55), (10, 7.5, 8.5), - (2, 1.5, 2.5)] + parameter_values = [ + (20, 15, 25), + (50, 45, 55), + (10, 7.5, 8.5), + (2, 1.5, 2.5), + ] # Fill parameter values and bounds parameters = [ - Mock(spec=bumps.parameter.Parameter, value=value, - bounds=Mock(limits=[lb, ub])) + Mock( + spec=bumps.parameter.Parameter, + value=value, + bounds=Mock(limits=[lb, ub]), + ) for value, lb, ub in parameter_values ] model = Mock(spec=refl1d.experiment.Experiment, xi=parameters) @@ -123,7 +137,7 @@ def test_fisher_workflow_refnx(refnx_model): [5.17704306e-06, 2.24179068e-06, -5.02221954e-07, -7.91886209e-07], [2.24179068e-06, 1.00559528e-06, -2.09433754e-07, -3.18583142e-07], [-5.02221954e-07, -2.09433754e-07, 5.75647233e-08, 1.03142100e-07], - [-7.91886209e-07, -3.18583142e-07, 1.03142100e-07, 1.99470835e-07] + [-7.91886209e-07, -3.18583142e-07, 1.03142100e-07, 1.99470835e-07], ] np.testing.assert_allclose(g, expected_fisher, rtol=1e-08) @@ -138,14 +152,15 @@ def test_fisher_workflow_refl1d(refl1d_model): [4.58294661e-06, 2.07712766e-06, -4.23068571e-07, -6.80596824e-07], [2.07712766e-06, 9.76175381e-07, -1.84017555e-07, -2.83513452e-07], [-4.23068571e-07, -1.84017555e-07, 4.51142562e-08, 8.21397190e-08], - [-6.80596824e-07, -2.83513452e-07, 8.21397190e-08, 1.62625881e-07] + [-6.80596824e-07, -2.83513452e-07, 8.21397190e-08, 1.62625881e-07], ] np.testing.assert_allclose(g, expected_fisher, rtol=1e-08) @patch('hogben.utils.reflectivity') -@pytest.mark.parametrize('model_class', ("mock_refl1d_model", - "mock_refnx_model")) +@pytest.mark.parametrize( + 'model_class', ('mock_refl1d_model', 'mock_refnx_model') +) def test_fisher_analytical_values(mock_reflectivity, model_class, request): """ Tests that the values of the calculated Fisher information matrix (FIM) @@ -197,15 +212,16 @@ def test_fisher_analytical_values(mock_reflectivity, model_class, request): g_correct = [ [1.28125, 0.5125, 25.625], [0.5125, 0.205, 10.25], - [25.625, 10.25, 512.5] + [25.625, 10.25, 512.5], ] g_reference = fisher(QS, xi, COUNTS, [model]) np.testing.assert_allclose(g_reference, g_correct, rtol=1e-08) @patch('hogben.utils.reflectivity') -@pytest.mark.parametrize('model_class', ("mock_refl1d_model", - "mock_refnx_model")) +@pytest.mark.parametrize( + 'model_class', ('mock_refl1d_model', 'mock_refnx_model') +) def test_fisher_importance_scaling(mock_reflectivity, model_class, request): """ Tests that the values of the calculated Fisher information matrix @@ -233,14 +249,13 @@ def test_fisher_importance_scaling(mock_reflectivity, model_class, request): g_correct = [ [1.28125, 1.025, 76.875], [0.5125, 0.41, 30.75], - [25.625, 20.5, 1537.5] + [25.625, 20.5, 1537.5], ] g_reference = fisher(QS, xi, COUNTS, [model]) np.testing.assert_allclose(g_reference, g_correct, rtol=1e-08) -@pytest.mark.parametrize('model_class', ("refl1d_model", - "refnx_model")) +@pytest.mark.parametrize('model_class', ('refl1d_model', 'refnx_model')) @pytest.mark.parametrize('step', (0.01, 0.0075, 0.0025, 0.001, 0.0001)) def test_fisher_consistent_steps(step, model_class, request): """ @@ -254,8 +269,9 @@ def test_fisher_consistent_steps(step, model_class, request): @patch('hogben.utils.reflectivity') -@pytest.mark.parametrize('model_class', ("mock_refl1d_model", - "mock_refnx_model")) +@pytest.mark.parametrize( + 'model_class', ('mock_refl1d_model', 'mock_refnx_model') +) @pytest.mark.parametrize('model_params', (1, 2, 3, 4)) def test_fisher_shape(mock_reflectivity, model_params, model_class, request): """ @@ -273,25 +289,33 @@ def test_fisher_shape(mock_reflectivity, model_params, model_class, request): @patch('hogben.utils.reflectivity') -@pytest.mark.parametrize('model_class', ("mock_refl1d_model", - "mock_refnx_model")) -@pytest.mark.parametrize('qs', - (np.arange(0.001, 1.0, 0.25), - np.arange(0.001, 1.0, 0.10), - np.arange(0.001, 1.0, 0.05), - np.arange(0.001, 1.0, 0.01))) -def test_fisher_diagonal_non_negative(mock_reflectivity, qs, model_class, - request): +@pytest.mark.parametrize( + 'model_class', ('mock_refl1d_model', 'mock_refnx_model') +) +@pytest.mark.parametrize( + 'qs', + ( + np.arange(0.001, 1.0, 0.25), + np.arange(0.001, 1.0, 0.10), + np.arange(0.001, 1.0, 0.05), + np.arange(0.001, 1.0, 0.01), + ), +) +def test_fisher_diagonal_non_negative( + mock_reflectivity, qs, model_class, request +): """Tests whether the diagonal values in the Fisher information matrix - are all zero or greater""" + are all zero or greater""" model = request.getfixturevalue(model_class) mock_reflectivity.side_effect = (np.random.rand(len(qs)) for _ in range(9)) counts = [np.ones(len(qs)) * 100] g = fisher([qs], model.xi, counts, [model]) assert np.all(np.diag(g)) >= 0 -@pytest.mark.parametrize('model_class', ("mock_refl1d_model", - "mock_refnx_model")) + +@pytest.mark.parametrize( + 'model_class', ('mock_refl1d_model', 'mock_refnx_model') +) @pytest.mark.parametrize('model_params', (1, 2, 3, 4)) def test_fisher_no_data(model_params, model_class, request): """Tests whether a model with zero data points properly returns an empty @@ -302,8 +326,9 @@ def test_fisher_no_data(model_params, model_class, request): np.testing.assert_equal(g, np.zeros((len(xi), len(xi)))) -@pytest.mark.parametrize('model_class', ("mock_refl1d_model", - "mock_refnx_model")) +@pytest.mark.parametrize( + 'model_class', ('mock_refl1d_model', 'mock_refnx_model') +) @patch('hogben.utils.reflectivity') def test_fisher_no_parameters(mock_reflectivity, model_class, request): """Tests whether a model without any parameters properly returns a @@ -314,8 +339,7 @@ def test_fisher_no_parameters(mock_reflectivity, model_class, request): np.testing.assert_equal(g.shape, (0, 0)) -@pytest.mark.parametrize('model_class', ("refnx_model", - "refl1d_model")) +@pytest.mark.parametrize('model_class', ('refnx_model', 'refl1d_model')) def test_fisher_doubling_with_two_identical_models(model_class, request): """ Tests that using two identical models with the same q-points and counts @@ -331,8 +355,7 @@ def test_fisher_doubling_with_two_identical_models(model_class, request): np.testing.assert_allclose(g_double, g_single * 2, rtol=1e-08) -@pytest.mark.parametrize('model_class', ("refnx_model", - "refl1d_model")) +@pytest.mark.parametrize('model_class', ('refnx_model', 'refl1d_model')) def test_multiple_models_shape(model_class, request): """ Tests that shape of the Fisher information matrix is equal to the total diff --git a/hogben/tests/test_simulation.py b/hogben/tests/test_simulation.py index 8f0e51dd..97a2a3e2 100644 --- a/hogben/tests/test_simulation.py +++ b/hogben/tests/test_simulation.py @@ -1,4 +1,3 @@ -import unittest import pytest import numpy as np @@ -29,14 +28,19 @@ class Test_Simulate(): model_1, data_1 = simulate(sample_1, angle_times, scale, bkg, dq, ref) def test_data_streaming(self): - """Tests that without an input for the datafile, the correct one is picked up""" + """Tests that without an input for the datafile, the correct one is + picked up""" angle_times = [(0.3, 100, 1000)] - _, simulated_datapoints = simulate(self.sample_1, angle_times, self.scale, self.bkg, self.dq, self.ref) - np.testing.assert_array_less(np.zeros(len(simulated_datapoints)), simulated_datapoints[:, 3]) # counts - - _, simulated_datapoints = simulate(self.sample_1, angle_times, self.scale, self.bkg, self.dq) - np.testing.assert_array_less(np.zeros(len(simulated_datapoints)), simulated_datapoints[:, 3]) # counts + _, simulated_datapoints = simulate( + self.sample_1, angle_times, self.scale, self.bkg, self.dq, + self.ref) + np.testing.assert_array_less(np.zeros(len(simulated_datapoints)), + simulated_datapoints[:, 3]) # counts + _, simulated_datapoints = simulate(self.sample_1, angle_times, + self.scale, self.bkg, self.dq) + np.testing.assert_array_less(np.zeros(len(simulated_datapoints)), + simulated_datapoints[:, 3]) # counts def test_refnx_simulate_model(self): """ @@ -58,11 +62,12 @@ def test_refnx_simulate_data(self): _, simulated_datapoints = simulate(self.sample_1, angle_times, self.scale, self.bkg, self.dq, self.ref) - + # reflectivity np.testing.assert_array_less(np.zeros(len(simulated_datapoints)), - simulated_datapoints[:,1]) # reflectivity + simulated_datapoints[:, 1]) + # counts np.testing.assert_array_less(np.zeros(len(simulated_datapoints)), - simulated_datapoints[:, 3]) # counts + simulated_datapoints[:, 3]) @pytest.mark.parametrize('instrument', ('OFFSPEC', @@ -84,18 +89,16 @@ def test_simulation_instruments(self, instrument): np.testing.assert_array_less(np.zeros(angle_times[0][1]), simulated_datapoints[:, 3]) # counts - @pytest.mark.parametrize('instrument', - ('OFFSPEC', - 'POLREF')) + @pytest.mark.parametrize('instrument', ('OFFSPEC', 'POLREF')) def test_simulation_magnetic_instruments(self, instrument): """ Tests that all of the instruments are able to simulate a model and counts data. """ angle_times = [(0.3, 100, 1000)] - _, simulated_datapoints = simulate_magnetic(self.sample_1, angle_times, - self.scale, self.bkg, self.dq, - inst_or_path=instrument) + _, simulated_datapoints = simulate_magnetic( + self.sample_1, angle_times, + self.scale, self.bkg, self.dq, inst_or_path=instrument) for i in range(4): # reflectivity diff --git a/hogben/underlayers.py b/hogben/underlayers.py index 81520083..daf06112 100644 --- a/hogben/underlayers.py +++ b/hogben/underlayers.py @@ -1,5 +1,4 @@ import os -import sys import time import numpy as np @@ -15,7 +14,7 @@ def _underlayer_results_visualise(save_path): save_path (str): path to directory to save results to. """ - from models.bilayers import BilayerDMPC, BilayerDPPC + from models.bilayers import BilayerDMPC # Choose sample here. bilayer = BilayerDMPC() @@ -38,8 +37,8 @@ def _underlayer_results_visualise(save_path): save_path, label) angle_times = [(0.7, 100, 40)] - underlayers = [(127.1, 5.39)] # Optimal DMPC bilayer underlayer. - #underlayers = [(76.5, 9.00)] # Optimal DPPC/Ra LPS bilayer underlayer. + underlayers = [(127.1, 5.39)] # Optimal DMPC bilayer underlayer. + # underlayers = [(76.5, 9.00)] # Optimal DPPC/Ra LPS bilayer underlayer. # Use nested sampling to validate the improvements. bilayer.nested_sampling([-0.56, 6.36], angle_times, save_path, @@ -48,6 +47,7 @@ def _underlayer_results_visualise(save_path): bilayer.nested_sampling([-0.56, 6.36], angle_times, save_path, 'D2O_H2O_with_underlayer', underlayers=underlayers) + def _underlayer_results_optimise(save_path): """Optimises choice of underlayer thickness and SLD for a bilayer sample. @@ -55,7 +55,7 @@ def _underlayer_results_optimise(save_path): save_path (str): path to directory to save results to. """ - from bilayers import BilayerDMPC, BilayerDPPC + from bilayers import BilayerDMPC # Choose sample here. bilayer = BilayerDMPC() @@ -74,7 +74,7 @@ def _underlayer_results_optimise(save_path): save_path = os.path.join(save_path, bilayer.name) file_path = os.path.join(save_path, 'optimised_underlayers.txt') with open(file_path, 'w') as file: - optimiser = Optimiser(bilayer) # Optimiser for the experiment. + optimiser = Optimiser(bilayer) # Optimiser for the experiment. # Calculate the objective value using no underlayers. g = optimiser.sample.underlayer_info(angle_times, contrasts, []) @@ -88,7 +88,8 @@ def _underlayer_results_optimise(save_path): # Calculate the objective value using a gold underlayer. underlayer = [(100, 4.7)] - g = optimiser.sample.underlayer_info(angle_times, contrasts, underlayer) + g = optimiser.sample.underlayer_info(angle_times, contrasts, + underlayer) val = -np.linalg.eigvalsh(g)[0] val = np.format_float_positional(val, precision=4, unique=False, fractional=False, trim='k') @@ -99,7 +100,8 @@ def _underlayer_results_optimise(save_path): # Calculate the objective value using a Permalloy underlayer. underlayer = [(100, 8.4)] - g = optimiser.sample.underlayer_info(angle_times, contrasts, underlayer) + g = optimiser.sample.underlayer_info(angle_times, contrasts, + underlayer) val = -np.linalg.eigvalsh(g)[0] val = np.format_float_positional(val, precision=4, unique=False, fractional=False, trim='k') @@ -127,11 +129,16 @@ def _underlayer_results_optimise(save_path): fractional=False, trim='k') # Save the conditions, objective value and computation time. - file.write('------------ {} Underlayers ------------\n'.format(num_underlayers)) - file.write('Thicknesses: {}\n'.format(list(np.round(thicknesses, 1)))) + file.write('------------ {} Underlayers ------------\n'.format( + num_underlayers)) + file.write('Thicknesses: {}\n'.format( + list(np.round(thicknesses, 1)))) + file.write('SLDs: {}\n'.format(list(np.round(slds, 2)))) file.write('Objective value: {}\n'.format(val)) - file.write('Computation time: {}\n\n'.format(round(end-start, 1))) + file.write('Computation time: {}\n\n'.format( + round(end - start, 1))) + if __name__ == '__main__': save_path = './results' diff --git a/hogben/utils.py b/hogben/utils.py index 952bf8e8..5d84eed3 100644 --- a/hogben/utils.py +++ b/hogben/utils.py @@ -223,7 +223,7 @@ def fisher(qs: list[np.ndarray], # scale with one if no importance was specified. importance_array = [] for param in xi: - if hasattr(param, "importance"): + if hasattr(param, 'importance'): importance_array.append(param.importance) else: importance_array.append(1) diff --git a/hogben/visualise.py b/hogben/visualise.py index 6f9bd34d..789f740f 100644 --- a/hogben/visualise.py +++ b/hogben/visualise.py @@ -1,5 +1,4 @@ import os -import sys import numpy as np import matplotlib.pyplot as plt @@ -7,14 +6,28 @@ from itertools import combinations from hogben.utils import save_plot -from hogben.models.base import VariableAngle, VariableContrast, VariableUnderlayer +from hogben.models.base import ( + BaseLipid, + BaseSample, + VariableAngle, + VariableContrast, + VariableUnderlayer, +) plt.rcParams['figure.figsize'] = (9, 7) plt.rcParams['figure.dpi'] = 600 -def angle_choice(sample, initial_angle_times, angle_range, points_new, - time_new, save_path, filename, contrasts=[]): +def angle_choice( + sample: BaseSample, + initial_angle_times: list, + angle_range: np.ndarray, + points_new: int, + time_new: type, + save_path: str, + filename: str, + contrasts: list = [], +): """Plots the minimum eigenvalue of the Fisher information matrix as a function of measurement angle. @@ -50,7 +63,7 @@ def angle_choice(sample, initial_angle_times, angle_range, points_new, g_new = sample.angle_info(new_angle_times, contrasts) # Combine the new information with the existing information. - min_eigs.append(np.linalg.eigvalsh(g_init+g_new)[0]) + min_eigs.append(np.linalg.eigvalsh(g_init + g_new)[0]) # Plot measurement angle versus minimum eigenvalue. fig = plt.figure() @@ -61,13 +74,22 @@ def angle_choice(sample, initial_angle_times, angle_range, points_new, # Save the plot. save_path = os.path.join(save_path, sample.name) - save_plot(fig, save_path, 'angle_choice_'+filename) + save_plot(fig, save_path, 'angle_choice_' + filename) # Return the angle with largest minimum eigenvalue. return angle_range[np.argmax(min_eigs)] -def angle_choice_with_time(sample, initial_angle, angle_range, time_range, - points, new_time, save_path, contrasts=[]): + +def angle_choice_with_time( + sample: BaseSample, + initial_angle: float, + angle_range: type, + time_range: type, + points: int, + new_time: float, + save_path: str, + contrasts: list = [], +): """Investigates how the second choice of angle for a `sample` changes as the counting time of the initial angle is increased. @@ -93,7 +115,7 @@ def angle_choice_with_time(sample, initial_angle, angle_range, time_range, ax.set_xlim(angle_range[0], angle_range[-1]) # Create the line that will have data added to it. - line, = ax.plot([], [], lw=3) + (line,) = ax.plot([], [], lw=3) # Initialiser function for the line. def init(): @@ -116,7 +138,7 @@ def animate(i): # Combine the information from the first and second angles. angle_times_new = [(new_angle, points, new_time)] g_new = sample.angle_info(angle_times_new, contrasts) - min_eigs.append(np.linalg.eigvalsh(g_init+g_new)[0]) + min_eigs.append(np.linalg.eigvalsh(g_init + g_new)[0]) # Update the data of the line. ax.set_ylim(min(min_eigs), max(min_eigs)) @@ -125,21 +147,28 @@ def animate(i): return line, # Define the animation. - anim_length = 4000 # Animation length in milliseconds. - frames = len(time_range) # Number of frames for animation. + anim_length = 4000 # Animation length in milliseconds. + frames = len(time_range) # Number of frames for animation. anim = FuncAnimation(fig, animate, init_func=init, blit=True, - frames=frames, interval=anim_length//frames) + frames=frames, interval=anim_length // frames) plt.close() # Save the animation as a .gif file. writergif = PillowWriter() save_path = os.path.join(save_path, sample.name, 'angle_choice_with_time.gif') + anim.save(save_path, writer=writergif) return anim -def contrast_choice_single(sample, contrast_range, initial_contrasts, - angle_times, save_path, filename): + +def contrast_choice_single(sample: BaseLipid, + contrast_range: np.ndarray, + initial_contrasts: list, + angle_times: list, + save_path: str, + filename: str + ) -> float: """Plots the minimum eigenvalue of the Fisher information matrix as a function of a single contrast SLD. @@ -170,7 +199,7 @@ def contrast_choice_single(sample, contrast_range, initial_contrasts, # Get the information from the new contrast and combine with initial. g_new = sample.contrast_info(angle_times, [new_contrast]) - min_eigs.append(np.linalg.eigvalsh(g_init+g_new)[0]) + min_eigs.append(np.linalg.eigvalsh(g_init + g_new)[0]) # Plot contrast SLD versus minimum eigenvalue. fig = plt.figure() @@ -182,11 +211,12 @@ def contrast_choice_single(sample, contrast_range, initial_contrasts, # Save the plot. save_path = os.path.join(save_path, sample.name) - save_plot(fig, save_path, 'contrast_choice_single_'+filename) + save_plot(fig, save_path, 'contrast_choice_single_' + filename) # Return the contrast SLD with largest minimum eigenvalue. return contrast_range[np.argmax(min_eigs)] + def contrast_choice_double(sample, contrast_range, angle_times, save_path): """Plots the minimum eigenvalue of the Fisher information matrix as a function of two contrast SLDs, assuming no prior measurement. @@ -221,12 +251,12 @@ def contrast_choice_double(sample, contrast_range, angle_times, save_path): # Duplicate the data so that the plot is not half-empty (or half-full?). # This is valid since the choice of contrast is commutative. # I.e., the order in which contrasts are measured is not important. - x = np.concatenate([contrasts[:,0], contrasts[:,1]]) - y = np.concatenate([contrasts[:,1], contrasts[:,0]]) + x = np.concatenate([contrasts[:, 0], contrasts[:, 1]]) + y = np.concatenate([contrasts[:, 1], contrasts[:, 0]]) min_eigs.extend(min_eigs) # Create a 3D plot of contrast pair versus minimum eigenvalue. - fig = plt.figure(figsize=[12,9]) + fig = plt.figure(figsize=[12, 9]) ax = fig.add_subplot(111, projection='3d') # Create the surface plot and add colour bar. @@ -240,7 +270,7 @@ def contrast_choice_double(sample, contrast_range, angle_times, save_path): ax.set_xlabel(x_label, fontsize=11, weight='bold') ax.set_ylabel(y_label, fontsize=11, weight='bold') ax.set_zlabel(z_label, fontsize=11, weight='bold') - ax.ticklabel_format(axis='z', style='sci', scilimits=(0,0)) + ax.ticklabel_format(axis='z', style='sci', scilimits=(0, 0)) # Save the plot. save_path = os.path.join(save_path, sample.name) @@ -261,8 +291,16 @@ def contrast_choice_double(sample, contrast_range, angle_times, save_path): maximum = np.argmax(min_eigs) return x[maximum], y[maximum] -def underlayer_choice(sample, thickness_range, sld_range, contrasts, - angle_times, save_path, filename=''): + +def underlayer_choice( + sample: BaseLipid, + thickness_range: np.ndarray, + sld_range: np.ndarray, + contrasts: list, + angle_times: list, + save_path: str, + filename: str = '', +): """Plots the minimum eigenvalue of the Fisher information matrix as a function of a sample's underlayer thickness and SLD. @@ -284,11 +322,11 @@ def underlayer_choice(sample, thickness_range, sld_range, contrasts, # Iterate over each underlayer thickness to investigate. x, y, min_eigs = [], [], [] - n = len(thickness_range)*len(sld_range) # Number of calculations. + n = len(thickness_range) * len(sld_range) # Number of calculations. for i, thick in enumerate(thickness_range): # Display progress. if i % 5 == 0: - print('>>> {0}/{1}'.format(i*len(sld_range), n)) + print('>>> {0}/{1}'.format(i * len(sld_range), n)) # Iterate over each underlayer SLD to investigate. for sld in sld_range: @@ -299,7 +337,7 @@ def underlayer_choice(sample, thickness_range, sld_range, contrasts, y.append(sld) # Create 3D plot of underlayer thickness and SLD versus minimum eigenvalue. - fig = plt.figure(figsize=[10,8]) + fig = plt.figure(figsize=[10, 8]) ax = fig.add_subplot(111, projection='3d') # Create the surface plot and add colour bar. @@ -313,11 +351,11 @@ def underlayer_choice(sample, thickness_range, sld_range, contrasts, ax.set_xlabel(x_label, fontsize=11, weight='bold') ax.set_ylabel(y_label, fontsize=11, weight='bold') ax.set_zlabel(z_label, fontsize=11, weight='bold') - ax.ticklabel_format(axis='z', style='sci', scilimits=(0,0)) + ax.ticklabel_format(axis='z', style='sci', scilimits=(0, 0)) # Save the plot. save_path = os.path.join(save_path, sample.name) - filename = 'underlayer_choice' + ('_'+filename if filename != '' else '') + filename = 'underlayer_choice' + ('_' + filename if filename != '' else '') save_plot(fig, save_path, filename) # Save different views of the 3D plot. diff --git a/requirements-with-tests.txt b/requirements-with-tests.txt index abc56a30..b4d0fe00 100644 --- a/requirements-with-tests.txt +++ b/requirements-with-tests.txt @@ -2,6 +2,7 @@ bumps corner dynesty emcee +flake8 importlib_resources matplotlib numpy diff --git a/setup.cfg b/setup.cfg index 5f1a8dc0..0d20c70d 100644 --- a/setup.cfg +++ b/setup.cfg @@ -33,3 +33,13 @@ install_requires = tqdm scipy include_package_data = True + +[flake8] +count = True +ignore = E741, W503, W605 +inline-quotes = single + +# `make_beam_spectra.py` uses a blank star import (F403), +# which makes it impossible to recognize certain modules (F405) +per-file-ignores = + make_beam_spectra.py: F403, F405