diff --git a/.github/workflows/publish.yaml b/.github/workflows/publish.yaml new file mode 100644 index 00000000..007b0cd5 --- /dev/null +++ b/.github/workflows/publish.yaml @@ -0,0 +1,206 @@ + +name: Publish + +on: + + push: + branches: [ master, main, develop, feature/github-actions] + tags: 'v*.*.*' + + pull_request: + branches: '*' + types: opened + + release: + types: [published] + + workflow_dispatch: + + repository_dispatch: + types: [release-event] + +jobs: + + init: + runs-on: ubuntu-latest + outputs: + build_type: ${{ steps.setvars.outputs.build_type }} + luna_version: ${{ steps.setversion.outputs.luna_version }} + + steps: + - name: Cancel previous workflow + uses: styfle/cancel-workflow-action@0.4.0 + with: + access_token: ${{ github.token }} + + - name: Set variables + id: setvars + #Only build the tests in a debug build, tests currently don't build in release mode + run: | + echo ${{github.base_ref}} + echo ${{github.ref}} + if [[ "${{github.base_ref}}" == "develop" || "${{github.ref}}" == "refs/heads/develop" ]]; then + echo "::set-output name=build_type::DEBUG" + echo Debug build + else + echo "::set-output name=build_type::RELEASE" + fi + + - name: Set Luna version + id: setversion + run: | + git clone https://github.com/project8/luna_base.git + cd luna_base + echo "::set-output name=luna_version::$(git describe --abbrev=0 --tags)" + + test_linux: + runs-on: ubuntu-latest + needs: [init] + + strategy: + matrix: + build_type: [DEBUG, RELEASE] + compiler: [gcc] + fail-fast: false + + steps: + - name: Checkout the repo + uses: actions/checkout@v2 + with: + submodules: recursive + + - name: Set env + run: | + if [[ "${{matrix.compiler}}" == "gcc" ]]; then + echo "CXX_VAL=g++" >> $GITHUB_ENV + else + echo "CXX_VAL=clang++" >> $GITHUB_ENV + fi + + - name: Print config + run: | + echo Build type: ${{ matrix.build_type }} + echo Compiler: ${{ matrix.compiler }} ${{ env.CXX_VAL }} + echo Luna version: ${{ needs.init.outputs.luna_version }} + + - name: Build with docker + run: | + cd ${{github.workspace}} + docker build \ + --build-arg IMG_USER=project8 \ + --build-arg IMG_REPO=p8compute_dependencies \ + --build-arg IMG_TAG=${{ needs.init.outputs.luna_version }} \ + --build-arg MERMITHID_TAG=test \ + --build-arg build_type=${{matrix.build_type}} \ + --build-arg CC_VAL=${{matrix.compiler}} \ + --build-arg CXX_VAL=${{env.CXX_VAL}} \ + --tag project8/mermithid:test \ + . + + publish: + if: github.event_name != 'pull_request' + name: Push Docker image to Docker Hub + runs-on: ubuntu-latest + needs: [init, test_linux] + steps: + + - name: Check out the repo + uses: actions/checkout@v2 + with: + submodules: recursive + + - name: Docker meta + id: docker_meta + uses: docker/metadata-action@v3 + with: + images: project8/mermithid + + - name: Set up QEMU + uses: docker/setup-qemu-action@v1 + + - name: Set up Docker Buildx + id: setup_buildx + uses: docker/setup-buildx-action@v1 + with: + buildkitd-flags: --debug + + - name: List platforms + run: | + echo Platforms: ${{ steps.setup_buildx.outputs.platforms }} + + - name: Log in to Docker Hub + uses: docker/login-action@v1 + with: + username: ${{ secrets.DOCKER_USERNAME }} + password: ${{ secrets.DOCKER_PASSWORD }} + + - name: Set env + run: | + echo "RELEASE_VERSION=${GITHUB_REF#refs/*/}" >> $GITHUB_ENV + + - name: Print config + run: | + echo Release version: ${{ env.RELEASE_VERSION }} + echo Luna version: ${{ needs.init.outputs.luna_version }} + + - name: Push to Docker Hub + id: build_push + uses: docker/build-push-action@v2 + with: + context: . + push: true + build-args: | + IMG_USER=project8 + IMG_REPO=p8compute_dependencies + IMG_TAG=${{ needs.init.outputs.luna_version }} + MERMITHID_TAG=${{ env.RELEASE_VERSION }} + build_type=${{needs.init.outputs.build_type}} + tags: ${{ steps.docker_meta.outputs.tags }} + +# disabled due to error +# doc: +# name: Build documentation +# runs-on: ubuntu-latest +# #needs: [init, test_linux] +# steps: +# +# - name: Check out the repo +# uses: actions/checkout@v2 +# with: +# submodules: recursive +# fetch-depth: 0 +# +# - name: Install dependencies +# run: | +# sudo apt-get update +# DEBIAN_FRONTEND=noninteractive sudo apt-get install -y \ +# tree \ +# doxygen \ +# python3-sphinx \ +# graphviz +# pip install sphinxcontrib-blockdiag +# pip install sphinxcontrib-contentui +# +# - name: Build docs +# run: | +# cd Documentation +# make + + notify-packages: + + if: (github.event_name == 'push' && contains(github.ref, 'refs/tags/') ) + || github.event_name == 'release' + || github.event_name == 'workflow_dispatch' + || github.event_name == 'repository_dispatch' + name: Trigger new docker P8Compute image + runs-on: ubuntu-latest + needs: [publish] + + steps: + + - name: Repository dispatch + uses: peter-evans/repository-dispatch@v1 + with: + token: ${{ secrets.REPO_ACCESS_TOKEN }} + repository: project8/luna + event-type: release-event diff --git a/Dockerfile b/Dockerfile index 91380fde..4922e950 100644 --- a/Dockerfile +++ b/Dockerfile @@ -7,9 +7,15 @@ FROM ${IMG_USER}/${IMG_REPO}:${IMG_TAG} as mermithid_common ARG build_type=Release ENV MERMITHID_BUILD_TYPE=$build_type -ENV MERMITHID_TAG=v1.2.2 +ARG MERMITHID_TAG=beta +ENV MERMITHID_TAG=${MERMITHID_TAG} ENV MERMITHID_BUILD_PREFIX=/usr/local/p8/mermithid/$MERMITHID_TAG +ARG CC_VAL=gcc +ENV CC=${CC_VAL} +ARG CXX_VAL=g++ +ENV CXX=${CXX_VAL} + RUN mkdir -p $MERMITHID_BUILD_PREFIX &&\ chmod -R 777 $MERMITHID_BUILD_PREFIX/.. &&\ cd $MERMITHID_BUILD_PREFIX &&\ @@ -22,6 +28,10 @@ RUN mkdir -p $MERMITHID_BUILD_PREFIX &&\ echo 'export PYTHONPATH=$MERMITHID_BUILD_PREFIX/$(python3 -m site --user-site | sed "s%$(python3 -m site --user-base)%%"):$PYTHONPATH' >> setup.sh &&\ /bin/true +RUN source $COMMON_BUILD_PREFIX/setup.sh &&\ + pip install iminuit &&\ + /bin/true + ######################## FROM mermithid_common as mermithid_done diff --git a/mermithid/misc/ComplexLineShapeUtilities.py b/mermithid/misc/ComplexLineShapeUtilities.py index 074013d0..4fe05dc7 100644 --- a/mermithid/misc/ComplexLineShapeUtilities.py +++ b/mermithid/misc/ComplexLineShapeUtilities.py @@ -58,6 +58,10 @@ def aseev_func_tail(energy_loss_array, gas_type): A2, omeg2, eps2 = 0.195, 14.13, 10.60 elif gas_type=="Kr": A2, omeg2, eps2 = 0.4019, 22.31, 16.725 + elif gas_type=="He": + A2, omeg2, eps2 = 0.1187, 33.40, 10.43 + elif gas_type=="Ar": + A2, omeg2, eps2 = 0.3344, 21.91, 21.14 return A2*omeg2**2./(omeg2**2.+4*(energy_loss_array-eps2)**2.) #convert oscillator strength into energy loss spectrum diff --git a/mermithid/misc/FakeTritiumDataFunctions.py b/mermithid/misc/FakeTritiumDataFunctions.py index bbf8edf2..5e36fe6e 100644 --- a/mermithid/misc/FakeTritiumDataFunctions.py +++ b/mermithid/misc/FakeTritiumDataFunctions.py @@ -7,6 +7,7 @@ from __future__ import absolute_import import numpy as np +import json from scipy.special import gamma from scipy import integrate @@ -18,7 +19,6 @@ from mermithid.misc.Constants import * from mermithid.misc.ConversionFunctions import * -from mermithid.misc.KrLineshapeFunctions import * """ Constants and functions used by processors/TritiumSpectrum/FakeDataGenerator.py @@ -147,13 +147,42 @@ def spectral_rate_in_window(K, Q, mnu, Kmin): else: return 0. + +def beta_rates(K, Q, mnu, index): + beta_rates = np.zeros(len(K)) + nu_mass_shape = ((Q - K[index])**2 -mnu**2)**0.5 + beta_rates[index] = GF**2.*Vud**2*Mnuc2/(2.*np.pi**3)*ephasespace(K[index], Q)*(Q - K[index])*nu_mass_shape + return beta_rates + + + # Unsmeared eta spectrum without a lower energy bound -def spectral_rate(K, Q, mnu): - if Q-mnu > K > 0: - return GF**2.*Vud**2*Mnuc2/(2.*np.pi**3)*ephasespace(K, Q)*(Q - K)*np.sqrt((Q - K)**2 - (mnu)**2) +def spectral_rate(K, Q, mnu, final_state_array): + + if isinstance(K, list) or isinstance(K, np.ndarray): + N_states = len(final_state_array[0]) + beta_rates_array = np.zeros([N_states, len(K)]) + + Q_states = Q+final_state_array[0]-np.max(final_state_array[0]) + + index = [np.where(K < Q_states[i]-mnu) for i in range(N_states)] + + beta_rates_array = [beta_rates(K, Q_states[i], mnu, index[i])*final_state_array[1][i] for i in range(N_states)] + to_return = np.nansum(beta_rates_array, axis=0)/np.nansum(final_state_array[1]) + + return to_return + else: - return 0. - #np.heaviside(Q-mnu-K, 0.5)*np.heaviside(K-V0, 0.5) + + return_value = 0. + + for i, e_binding in enumerate(final_state_array[0]): + # binding energies are negative + Q_state = Q+e_binding + if Q_state-mnu > K > 0: + return_value += final_state_array[1][i] *(GF**2.*Vud**2*Mnuc2/(2.*np.pi**3)*ephasespace(K, Q_state)*(Q_state - K)*np.sqrt((Q_state - K)**2 - (mnu)**2)) + + return return_value/np.sum(final_state_array[1]) #Flat background with lower and upper bounds Kmin and Kmax @@ -169,12 +198,12 @@ def bkgd_rate(): return 1. -##Lineshape option: In this case, simply a Gaussian -#def gaussian(x,a): -# return 1/((2.*np.pi)**0.5*a[0])*(np.exp(-0.5*((x-a[1])/a[0])**2)) +#Lineshape option: In this case, simply a Gaussian +def gaussian(x,a): + return 1/((2.*np.pi)**0.5*a[0])*(np.exp(-0.5*((x-a[1])/a[0])**2)) -#Normalized simplified linesahape with scattering +#Normalized simplified lineshape with scattering def simplified_ls(K, Kcenter, FWHM, prob, p0, p1, p2, p3): sig0 = FWHM/float(2*np.sqrt(2*np.log(2))) shape = gaussian(K, [sig0, Kcenter]) @@ -245,7 +274,8 @@ def convolved_bkgd_rate(K, Kmin, Kmax, lineshape, ls_params, min_energy, max_ene #Convolution of signal and lineshape using scipy.signal.convolve def convolved_spectral_rate_arrays(K, Q, mnu, Kmin, - lineshape, ls_params, min_energy, max_energy): + lineshape, ls_params, min_energy, max_energy, + complexLineShape, final_state_array): """K is an array-like object """ logger.info('Using scipy convolve') @@ -263,11 +293,12 @@ def convolved_spectral_rate_arrays(K, Q, mnu, Kmin, elif lineshape=='simplified_scattering' or lineshape=='simplified': lineshape_rates = simplified_ls(K_lineshape, 0, ls_params[0], ls_params[1], ls_params[2], ls_params[3], ls_params[4], ls_params[5]) elif lineshape=='detailed_scattering' or lineshape=='detailed': - lineshape_rates = spectrum_func(K_lineshape/1000., ls_params[0], 0, ls_params[1], 1) - beta_rates = np.zeros(len(K)) - for i,ke in enumerate(K): - beta_rates[i] = spectral_rate(ke, Q, mnu) + lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, 1, ls_params[1]) + + beta_rates = spectral_rate(K, Q, mnu, final_state_array) #np.zeros(len(K)) + #for i,ke in enumerate(K): + # beta_rates[i] = spectral_rate(ke, Q, mnu, final_state_array) #Convolving convolved = convolve(beta_rates, lineshape_rates, mode='same') @@ -278,7 +309,7 @@ def convolved_spectral_rate_arrays(K, Q, mnu, Kmin, #Convolution of background and lineshape using scipy.signal.convolve -def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy, max_energy): +def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy, max_energy, complexLineShape): """K is an array-like object """ energy_half_range = max(max_energy, abs(min_energy)) @@ -294,7 +325,7 @@ def convolved_bkgd_rate_arrays(K, Kmin, Kmax, lineshape, ls_params, min_energy, elif lineshape=='simplified_scattering' or lineshape=='simplified': lineshape_rates = simplified_ls(K_lineshape, 0, ls_params[0], ls_params[1], ls_params[2], ls_params[3], ls_params[4], ls_params[5]) elif lineshape=='detailed_scattering' or lineshape=='detailed': - lineshape_rates = spectrum_func(K_lineshape/1000., ls_params[0], 0, ls_params[1], 1) + lineshape_rates = complexLineShape.spectrum_func_1(K_lineshape/1000., ls_params[0], 0, 1, ls_params[1]) bkgd_rates = np.full(len(K), bkgd_rate()) if len(K) < len(K_lineshape): diff --git a/mermithid/misc/KrLineshapeFunctions.py b/mermithid/misc/KrLineshapeFunctions.py deleted file mode 100644 index 4fa791dd..00000000 --- a/mermithid/misc/KrLineshapeFunctions.py +++ /dev/null @@ -1,318 +0,0 @@ -''' -Various functions for the Krypton lineshape computation -Adapted from https://github.com/project8/scripts/blob/feature/KrScatterFit/machado/fitting/data_fitting_rebuild.py#L429 - by E. M. Machado -Changes from that version: -- Kr specific functions (shake-up/shake-off) removed -- Spectrum is normalized - -Author: T. Weiss, C. Claessens, Y. Sun -Date:4/6/2020 -''' - -from __future__ import absolute_import - -import numpy as np - -from scipy import integrate -from scipy.signal import convolve -import os - -from morpho.utilities import morphologging -logger = morphologging.getLogger(__name__) - -from mermithid.misc.Constants import * -from mermithid.misc.ConversionFunctions import * - - -# Natural constants -kr_line = kr_k_line_e()*1e-3 #17.8260 # keV -kr_line_width = kr_k_line_width() #2.83 # eV -e_charge = e() #1.60217662*10**(-19) # Coulombs , charge of electron -m_e = m_electron()/e()*c()**2 #9.10938356*10**(-31) # Kilograms , mass of electron -mass_energy_electron = m_electron()*1e-3 #510.9989461 # keV -max_scatters = 20 -gases = ["H2"] - -# This is an important parameter which determines how finely resolved -# the scatter calculations are. 10000 seems to produce a stable fit, with minimal slowdown -num_points_in_std_array = 10000 - -# Establishes a standard energy loss array (SELA) from -1000 eV to 1000 eV -# with number of points equal to num_points_in_std_array. All convolutions -# will be carried out on this particular discretization -def std_eV_array(): - emin = -1000 - emax = 1000 - array = np.linspace(emin,emax,num_points_in_std_array) - return array - -# A lorentzian function -def lorentzian(x_array,x0,FWHM): - HWHM = FWHM/2. - func = (1./np.pi)*HWHM/((x_array-x0)**2.+HWHM**2.) - return func - -# A lorentzian line centered at 0 eV, with 2.83 eV width on the SELA -def std_lorentzian_17keV(): - x_array = std_eV_array() - ans = lorentzian(x_array,0,kr_line_width) - return ans - - -# A gaussian function -def gaussian(x_array,A,sigma=0,mu=0): - if isinstance(A, list): - a = A - x = x_array - return 1/((2.*np.pi)**0.5*a[0])*(np.exp(-0.5*((x-a[1])/a[0])**2)) - else: - f = A*(1./(sigma*np.sqrt(2*np.pi)))*np.exp(-(((x_array-mu)/sigma)**2.)/2.) - return f - -# A gaussian centered at 0 eV with variable width, on the SELA -def std_gaussian(sigma): - x_array = std_eV_array() - ans = gaussian(x_array,1,sigma,0) - return ans - -def std_dirac(): - x_array = std_eV_array() - ans = np.zeros(len(x_array)) - min_x = np.min(np.abs(x_array)) - ans[np.abs(x_array)==min_x] = 1. - logger.warning('Spectrum will be shifted by lineshape by {} eV'.format(min_x)) - if min_x > 0.1: - logger.warning('Lineshape will shift spectrum by > 0.1 eV') - if min_x > 1.: - logger.warning('Lineshape will shift spectrum by > 1 eV') - raise ValueError('problem with std_eV_array()') - return ans - -# Converts a gaussian's FWHM to sigma -def gaussian_FWHM_to_sigma(FWHM): - sigma = FWHM/(2*np.sqrt(2*np.log(2))) - return sigma - -# Converts a gaussian's sigma to FWHM -def gaussian_sigma_to_FWHM(sigma): - FWHM = sigma*(2*np.sqrt(2*np.log(2))) - return FWHM - -# normalizes a function, but depends on binning. -# Only to be used for functions evaluated on the SELA -def normalize(f): - x_arr = std_eV_array() - f_norm = integrate.simps(f,x=x_arr) - # if f_norm < 0.99 or f_norm > 1.01: - # print(f_norm) - f_normed = f/f_norm - return f_normed - -#returns array with energy loss/ oscillator strength data -#def read_oscillator_str_file(filename): -# f = open(filename, "r") -# lines = f.readlines() -# energyOsc = [[],[]] #2d array of energy losses, oscillator strengths -# -# for line in lines: -# if line != "" and line[0]!="#": -# raw_data = [float(i) for i in line.split("\t")] -# energyOsc[0].append(raw_data[0]) -# energyOsc[1].append(raw_data[1]) -# -# energyOsc = np.array(energyOsc) -# ### take data and sort by energy -# sorted_indices = energyOsc[0].argsort() -# energyOsc[0] = energyOsc[0][sorted_indices] -# energyOsc[1] = energyOsc[1][sorted_indices] -# return energyOsc - -# A sub function for the scatter function. Found in -# "Energy loss of 18 keV electrons in gaseous T and quench condensed D films" -# by V.N. Aseev et al. 2000 -#def aseev_func_tail(energy_loss_array, gas_type): -# if gas_type=="H2": -# A2, omeg2, eps2 = 0.195, 14.13, 10.60 -# elif gas_type=="Kr": -# A2, omeg2, eps2 = 0.4019, 22.31, 16.725 -# return A2*omeg2**2./(omeg2**2.+4*(energy_loss_array-eps2)**2.) - -#convert oscillator strength into energy loss spectrum -#def get_eloss_spec(e_loss, oscillator_strength, kinetic_en = kr_line * 1000): #energies in eV -# e_rydberg = 13.605693009 #rydberg energy (eV) -# a0 = 5.291772e-11 #bohr radius -# return np.where(e_loss>0 , 4.*np.pi*a0**2 * e_rydberg / (kinetic_en * e_loss) * oscillator_strength * np.log(4. * kinetic_en * e_loss / (e_rydberg**3.) ), 0) - -# Function for energy loss from a single scatter of electrons by -# V.N. Aseev et al. 2000 -# This function does the work of combining fit_func1 and fit_func2 by -# finding the point where they intersect. -# Evaluated on the SELA -#def single_scatter_f(gas_type): -# energy_loss_array = std_eV_array() -# f = 0 * energy_loss_array -# -# input_filename = "../data/" + gas_type + "OscillatorStrength.txt" -# energy_fOsc = read_oscillator_str_file(input_filename) -# fData = interpolate.interp1d(energy_fOsc[0], energy_fOsc[1], kind='linear') -# for i in range(len(energy_loss_array)): -# if energy_loss_array[i] < energy_fOsc[0][0]: -# f[i] = 0 -# elif energy_loss_array[i] <= energy_fOsc[0][-1]: -# f[i] = fData(energy_loss_array[i]) -# else: -# f[i] = aseev_func_tail(energy_loss_array[i], gas_type) -# -# f_e_loss = get_eloss_spec(energy_loss_array, f) -# f_normed = normalize(f_e_loss) -# #plt.plot(energy_loss_array, f_e_loss) -# #plt.show() -# return f_normed - -# Convolves a function with the single scatter function, on the SELA -#def another_scatter(input_spectrum, gas_type): -# single = single_scatter_f(gas_type) -# f = convolve(single,input_spectrum,mode='same') -# f_normed = normalize(f) -# return f_normed - -# Convolves the scatter function with itself over and over again and saves -# the results to .npy files. -#def generate_scatter_convolution_files(gas_type): -# t = time.time() -# first_scatter = single_scatter_f(gas_type) -# scatter_num_array = range(2,max_scatters+1) -# current_scatter = first_scatter -# np.save('scatter_spectra_files/scatter'+gas_type+"_"+str(1).zfill(2),current_scatter) -# # x = std_eV_array() # diagnostic -# for i in scatter_num_array: -# current_scatter = another_scatter(current_scatter, gas_type) -# np.save('scatter_spectra_files/scatter'+gas_type+"_"+str(i).zfill(2),current_scatter) -# # plt.plot(x,current_scatter) # diagnostic -# # plt.show() # diagnostic -# elapsed = time.time() - t -# logger.info('Files generated in '+str(elapsed)+'s') -# return - -# Returns the name of the current path -def get_current_path(): - path = os.path.abspath(os.getcwd()) - return path - -# Prints a list of the contents of a directory -def list_files(path): - directory = os.popen("ls "+path).readlines() - strippeddirs = [s.strip('\n') for s in directory] - return strippeddirs - -# Returns the name of the current directory -def get_current_dir(): - current_path = os.popen("pwd").readlines() - stripped_path = [s.strip('\n') for s in current_path] - stripped_path = stripped_path[0].split('/') - current_dir = stripped_path[len(stripped_path)-1] - return current_dir - -# Checks for the existence of a directory called 'scatter_spectra_files' -# and checks that this directory contains the scatter spectra files. -# If not, this function calls generate_scatter_convolution_files. -# This function also checks to make sure that the scatter files have the correct -# number of points in the SELA, and if not, it generates fresh files -#def check_existence_of_scatter_files(gas_type): -# current_path = get_current_path() -# current_dir = get_current_dir() -# stuff_in_dir = list_files(current_path) -# if 'scatter_spectra_files' not in stuff_in_dir and current_dir != 'scatter_spectra_files': -# logger.warning('Scatter files not found, generating') -# os.popen("mkdir scatter_spectra_files") -# time.sleep(2) -# generate_scatter_convolution_files(gas_type) -# else: -# directory = os.popen("ls scatter_spectra_files").readlines() -# strippeddirs = [s.strip('\n') for s in directory] -# if len(directory) != len(gases) * max_scatters: -# generate_scatter_convolution_files(gas_type) -# test_file = 'scatter_spectra_files/scatter'+gas_type+'_01.npy' -# test_arr = np.load(test_file) -# if len(test_arr) != num_points_in_std_array: -# logger.warning('Scatter files do not match standard array binning, generating fresh files') -# generate_scatter_convolution_files(gas_type) -# return - -# Given a function evaluated on the SELA, convolves it with a gaussian -def convolve_gaussian(func_to_convolve,gauss_FWHM_eV): - sigma = gaussian_FWHM_to_sigma(gauss_FWHM_eV) - resolution_f = std_gaussian(sigma) - ans = convolve(resolution_f,func_to_convolve,mode='same') - ans_normed = normalize(ans) - return ans_normed - -# Produces a full spectral shape on the SELA, given a gaussian resolution -# and a scatter probability -def make_spectrum(gauss_FWHM_eV,scatter_prob,gas_type, emitted_peak='dirac'): - current_path = get_current_path() - # check_existence_of_scatter_files() - #filenames = list_files('scatter_spectra_files') - scatter_num_array = range(1,max_scatters+1) - en_array = std_eV_array() - current_full_spectrum = np.zeros(len(en_array)) - if emitted_peak == 'lorentzian': - current_working_spectrum = std_lorentzian_17keV() #Normalized - elif emitted_peak == 'shake': - current_working_spectrum = shake_spectrum() - else: - current_working_spectrum = std_dirac() - current_working_spectrum = convolve_gaussian(current_working_spectrum,gauss_FWHM_eV) #Still normalized - zeroth_order_peak = current_working_spectrum - current_full_spectrum += current_working_spectrum #I believe, still normalized - norm = 1 - for i in scatter_num_array: - current_working_spectrum = np.load(os.path.join(current_path, 'scatter_spectra_files/scatter'+gas_type+'_'+str(i).zfill(2)+'.npy')) - current_working_spectrum = normalize(convolve(zeroth_order_peak,current_working_spectrum,mode='same')) - current_full_spectrum += current_working_spectrum*scatter_prob**scatter_num_array[i-1] - norm += scatter_prob**scatter_num_array[i-1] - # plt.plot(en_array,current_full_spectrum) # diagnostic - # plt.show() # diagnostic - return current_full_spectrum/norm - -# Takes only the nonzero bins of a histogram -def get_only_nonzero_bins(bins,hist): - nonzero_idx = np.where(hist!=0) - hist_nonzero = hist[nonzero_idx] - hist_err = np.sqrt(hist_nonzero) - bins_nonzero = bins[nonzero_idx] - return bins_nonzero , hist_nonzero , hist_err - -# Flips an array left-to-right. Useful for converting between energy and frequency -def flip_array(array): - flipped = np.fliplr([array]).flatten() - return flipped - - -def spectrum_func(x_keV, *p0): - logger.info('Using detailed lineshape') - x_eV = x_keV*1000. - en_loss_array = std_eV_array() - en_loss_array_min = en_loss_array[0] - en_loss_array_max = en_loss_array[len(en_loss_array)-1] - en_array_rev = flip_array(-1*en_loss_array) - f = np.zeros(len(x_keV)) - - FWHM_G_eV = p0[0] - line_pos_keV = p0[1] - - line_pos_eV = line_pos_keV*1000. - x_eV_minus_line = x_eV - line_pos_eV - zero_idx = np.r_[np.where(x_eV_minus_line<-1*en_loss_array_max)[0],np.where(x_eV_minus_line>-1*en_loss_array_min)[0]] - nonzero_idx = [i for i,_ in enumerate(x_keV) if i not in zero_idx] - - for gas_index in range(len(gases)): - gas_type = gases[gas_index] - scatter_prob = p0[2*gas_index+2] - amplitude = p0[2*gas_index+3] - - full_spectrum = make_spectrum(FWHM_G_eV,scatter_prob, gas_type) - full_spectrum_rev = flip_array(full_spectrum) - f[nonzero_idx] += amplitude*np.interp(x_eV_minus_line[nonzero_idx],en_array_rev,full_spectrum_rev) - return f - diff --git a/mermithid/misc/__init__.py b/mermithid/misc/__init__.py index 0dc77499..afb4ebbe 100644 --- a/mermithid/misc/__init__.py +++ b/mermithid/misc/__init__.py @@ -8,4 +8,3 @@ from . import TritiumFormFactor from . import FakeTritiumDataFunctions from . import ConversionFunctions -from . import KrLineshapeFunctions diff --git a/mermithid/misc/saenz_mfs.json b/mermithid/misc/saenz_mfs.json new file mode 100644 index 00000000..3d3b9d7d --- /dev/null +++ b/mermithid/misc/saenz_mfs.json @@ -0,0 +1 @@ +{"Binding energy": [1.897, 1.844, 1.773, 1.65, 1.546, 1.455, 1.341, 1.232, 1.138, 1.047, 0.96, 0.849, 0.754, 0.647999999999999, 0.538, 0.446, 0.345, 0.24, 0.151999999999999, 0.0629999999999999, -0.0429999999999999, -0.147, -0.247, -0.347, -0.446999999999999, -0.613, -0.865, -1.112, -1.36, -1.61, -1.86, -2.186, -2.68199999999999, -3.23499999999999, -3.75, -16.603, -17.603, -18.799, -19.761, -20.73, -21.701, -22.676, -23.653, -24.632, -25.613, -26.596, -27.581, -28.567, -29.558, -30.593, -31.66, -32.637, -33.595, -34.562, -35.548, -36.566, -37.602, -38.609, -39.601, -40.601, -41.607, -42.614, -43.597, -44.584, -45.586, -46.616, -47.601, -48.565, -49.604, -50.599, -51.594, -52.605, -53.611, -54.629, -55.621, -56.632, -57.621, -58.608, -59.608, -60.604, -61.133, -62.615, -63.607, -64.613, -65.604, -66.595, -67.592, -68.589, -69.5759999999999, -70.5809999999999, -71.589, -72.5459999999999, -73.5489999999999, -74.568, -75.533, -76.6149999999999, -77.5669999999999, -78.607, -79.613, -80.6259999999999, -81.6079999999999, -82.6019999999999, -83.5929999999999, -84.5939999999999, -85.601, -86.601, -87.598, -89.0699999999999, -91.086, -93.0849999999999, -95.0849999999999, -97.084, -99.084, -101.086, -103.087, -105.088, -107.089, -109.089, -111.09, -113.09, -115.09, -117.091, -119.091, -121.091, -123.092, -125.092, -127.092, -129.092, -131.093, -133.093, -135.093, -137.093, -140.065, -144.067, -148.068, -152.069, -156.07, -160.072, -164.073, -168.074, -172.075, -176.076, -180.076, -184.077, -188.078, -192.079, -196.079, -200.08, -204.08, -208.081, -212.082, -216.082, -220.082, -224.083, -228.083, -232.084, -236.084], "Probability": [9.99999999999999e-08, 0.0006900000000000001, 0.00046, 0.00233, 0.00553, 0.00457, 0.02033, 0.01649, 0.03877, 0.038079999999999996, 0.06809, 0.11214, 0.10112, 0.24406, 0.32337000000000005, 0.40864, 0.68745, 0.66279, 0.51412, 0.6556100000000001, 0.54588, 0.37231000000000003, 0.25473, 0.16959, 0.11369, 0.16946999999999998, 0.10094, 0.05732, 0.02806, 0.013160000000000002, 0.00623, 0.0042, 0.0008, 0.00015, 0.0, 0.0, 0.0, 1.1999999999999999e-05, 0.000113, 0.0006560000000000001, 0.002567, 0.007149, 0.014804, 0.023583, 0.029715, 0.030307, 0.025527, 0.018080000000000002, 0.01107, 0.007377000000000001, 0.010637, 0.019095, 0.022178, 0.016434, 0.009037, 0.004989, 0.003978, 0.004124, 0.004152, 0.0039250000000000005, 0.003457, 0.003186, 0.0027010000000000003, 0.0027129999999999997, 0.002481, 0.002412, 0.001907, 0.001938, 0.0017599999999999998, 0.001575, 0.0015409999999999998, 0.001485, 0.001557, 0.001895, 0.002427, 0.003357, 0.004095, 0.004714, 0.005033999999999999, 0.005152, 0.005442000000000001, 0.005859, 0.006617, 0.0070940000000000005, 0.007404, 0.007164, 0.006563, 0.005620000000000001, 0.004691, 0.00368, 0.003049, 0.00221, 0.001928, 0.0017610000000000002, 0.0015300000000000001, 0.001215, 0.0013900000000000002, 0.001216, 0.0014219999999999999, 0.001384, 0.001368, 0.001316, 0.001153, 0.0010760000000000001, 0.000921, 0.0007570000000000001, 0.000696, 0.0006180000000000001, 0.00054, 0.00048300000000000003, 0.00043200000000000004, 0.000388, 0.00035150000000000003, 0.00031800000000000003, 0.000289, 0.000264, 0.00024150000000000002, 0.000222, 0.0002045, 0.000189, 0.00017500000000000003, 0.00016250000000000002, 0.000151, 0.000141, 0.0001315, 0.000123, 0.000115, 0.00010800000000000001, 0.0001015, 9.549999999999999e-05, 8.999999999999999e-05, 8.45e-05, 7.775e-05, 6.95e-05, 6.25e-05, 5.625e-05, 5.1000000000000006e-05, 4.625e-05, 4.225e-05, 3.85e-05, 3.525e-05, 3.25e-05, 3e-05, 2.775e-05, 2.5750000000000002e-05, 2.375e-05, 2.225e-05, 2.075e-05, 1.925e-05, 1.8e-05, 1.7e-05, 1.6000000000000003e-05, 1.5e-05, 1.4e-05, 1.325e-05, 1.1750000000000001e-05, NaN]} \ No newline at end of file diff --git a/mermithid/processors/Fitters/BinnedDataFitter.py b/mermithid/processors/Fitters/BinnedDataFitter.py new file mode 100644 index 00000000..edd2faac --- /dev/null +++ b/mermithid/processors/Fitters/BinnedDataFitter.py @@ -0,0 +1,193 @@ +''' +Author: C. Claessens +Date:6/12/2020 +Description: + Minimizes poisson loglikelihood using iMinuit. + Configure with arbitrary histogram and model. + The model should consist of a pdf multiplied with a free amplitude parameter representing the total number of events. + +''' + +from __future__ import absolute_import + +import numpy as np +import scipy +import sys +import json +from iminuit import Minuit +from morpho.utilities import morphologging, reader +from morpho.processors import BaseProcessor + +logger = morphologging.getLogger(__name__) + + + +__all__ = [] +__all__.append(__name__) + +class BinnedDataFitter(BaseProcessor): + ''' + Processor that + Args: + variables: dictionary key under which data is stored + parameter_names: names of model parameters + initial_values: list of parameter initial values + limits: parameter limits given as [lower, upper] for each parameter + fixed: boolean list of same length as parameters. Only un-fixed parameters will be fitted + bins: bins for data histogram + binned_data: if True data is assumed to already be histogrammed + print_level: if 1 fit plan and result summaries are printed + constrained_parameter_indices: list of indices indicating which parameters are contraint. Contstraints will be Gaussian. + constrained_parameter_means: Mean of Gaussian constraint + constrained_parameter_widths: Standard deviation of Gaussian constrained + + Inputs: + data: + Output: + result: dictionary containing fit results and uncertainties + ''' + def InternalConfigure(self, params): + ''' + Configure + ''' + + # Read other parameters + self.namedata = reader.read_param(params, 'variables', "required") + self.parameter_names = reader.read_param(params, 'parameter_names', ['A', 'mu', 'sigma']) + self.initial_values = reader.read_param(params, 'initial_values', [1]*len(self.parameter_names)) + self.limits = reader.read_param(params, 'limits', [[None, None]]*len(self.parameter_names)) + self.fixes = reader.read_param(params, 'fixed', [False]*len(self.parameter_names)) + self.bins = reader.read_param(params, 'bins', np.linspace(-2, 2, 100)) + self.binned_data = reader.read_param(params, 'binned_data', False) + self.print_level = reader.read_param(params, 'print_level', 1) + self.constrained_parameters = reader.read_param(params, 'constrained_parameter_indices', []) + self.constrained_means = reader.read_param(params, 'constrained_parameter_means', []) + self.constrained_widths = reader.read_param(params, 'constrained_parameter_widths', []) + + # derived configurations + self.bin_centers = self.bins[0:-1]+0.5*(self.bins[1]-self.bins[0]) + self.parameter_errors = [max([0.1, 0.1*p]) for p in self.initial_values] + + return True + + def InternalRun(self): + logger.info('namedata: {}'.format(self.namedata)) + + if self.binned_data: + logger.info('Data is already binned. Getting binned data') + self.hist = self.data[self.namedata] + if len(self.hist) != len(self.bin_centers): + logger.error('Number of bins and histogram entries do not match') + raise ValueError('Number of bins and histogram entries do not match') + else: + logger.info('Data is unbinned. Will histogram before fitting.') + self.hist, _ = np.histogram(self.data[self.namedata], self.bins) + + result_array, error_array = self.fit() + + # save results + self.results = {} + self.results['param_values'] = result_array + self.results['param_errors'] = error_array + self.results['correlation_matrix'] = np.array(self.m_binned.covariance.correlation()) + for i, k in enumerate(self.parameter_names): + self.results[k] = {'value': result_array[i], 'error': error_array[i]} + + return True + + + + def model(self, x, A, mu, sigma): + """ + Overwrite by whatever Model + """ + f = A*(1/(sigma*np.sqrt(2*np.pi)))*np.exp(-(((x-mu)/sigma)**2.)/2.) + return f + + + + def fit(self): + # Now minimize neg log likelihood using iMinuit + if self.print_level > 0: + logger.info('This is the plan:') + logger.info('Fitting data consisting of {} elements'.format(np.sum(self.hist))) + logger.info('Fit parameters: {}'.format(self.parameter_names)) + logger.info('Initial values: {}'.format(self.initial_values)) + logger.info('Initial error: {}'.format(self.parameter_errors)) + logger.info('Limits: {}'.format(self.limits)) + logger.info('Fixed in fit: {}'.format(self.fixes)) + logger.info('Constrained parameters: {}'.format([self.parameter_names[i] for i in self.constrained_parameters])) + + m_binned = Minuit(self.negPoissonLogLikelihood, + self.initial_values, + # error=self.parameter_errors, + # errordef = 0.5, limit = self.limits, + name=self.parameter_names, + # fix=self.fixes, + # print_level=self.print_level, + # throw_nan=True + ) + + m_binned.errordef = 0.5 + m_binned.errors = self.parameter_errors + m_binned.throw_nan = True + m_binned.print_level = self.print_level + + for i, name in enumerate(self.parameter_names): + m_binned.fixed[name] = self.fixes[i] + m_binned.limits[name] = self.limits[i] + + # minimze + m_binned.simplex().migrad() + m_binned.hesse() + #self.param_states = m_binned.get_param_states() + self.m_binned = m_binned + + # results + result_array = np.array(m_binned.values) + error_array = np.array(m_binned.errors) + + if self.print_level == 1: + logger.info('Fit results: {}'.format(result_array)) + logger.info('Errors: {}'.format(error_array)) + + return result_array, error_array + + + def PoissonLogLikelihood(self, params): + + + # expectation + model_return = self.model(self.bin_centers, *params) + if np.shape(model_return)[0] == 2: + expectation, expectation_error = model_return + else: + expectation = model_return + + if np.min(expectation) < 0: + logger.error('Expectation contains negative numbers: Minimum {} -> {}. They will be excluded but something could be horribly wrong.'.format(np.argmin(expectation), np.min(expectation))) + logger.error('FYI, the parameters are: {}'.format(params)) + + # exclude bins where expectation is <= zero or nan + index = np.where(expectation>0) + + # poisson log likelihoood + ll = (self.hist[index]*np.log(expectation[index]) - expectation[index]).sum() + + # extended ll: poisson total number of events + N = np.nansum(expectation) + extended_ll = -N+np.sum(self.hist)*np.log(N)+ll + return extended_ll + + + def negPoissonLogLikelihood(self, params): + + # neg log likelihood + neg_ll = - self.PoissonLogLikelihood(params) + + # constrained parameters + if len(self.constrained_parameters) > 0: + for i, param in enumerate(self.constrained_parameters): + neg_ll += 0.5 * ((params[param] - self.constrained_means[i])/ self.constrained_widths[i])**2 + + return neg_ll diff --git a/mermithid/processors/Fitters/__init__.py b/mermithid/processors/Fitters/__init__.py new file mode 100644 index 00000000..634545ab --- /dev/null +++ b/mermithid/processors/Fitters/__init__.py @@ -0,0 +1,7 @@ +''' +''' + +from __future__ import absolute_import + +from .BinnedDataFitter import BinnedDataFitter + diff --git a/mermithid/processors/IO/MultiChannelCicadaReader.py b/mermithid/processors/IO/MultiChannelCicadaReader.py index e3b3bc3c..69241b17 100644 --- a/mermithid/processors/IO/MultiChannelCicadaReader.py +++ b/mermithid/processors/IO/MultiChannelCicadaReader.py @@ -99,10 +99,10 @@ def Reader(self): index = np.where((true_frequencies>=self.transition_freqs[i][0]) &(true_frequencies K conversion - #Simplified scattering model parameters + #Scattering model parameters self.survival_prob = reader.read_param(params, 'survival_prob', 0.77) self.scattering_sigma = reader.read_param(params, 'scattering_sigma', 18.6) self.NScatters = reader.read_param(params, 'NScatters', 20) - + self.scatter_proportion = reader.read_param(params, 'scatter_proportion', 1.0) + self.min_energy = reader.read_param(params,'min_lineshape_energy', -1000) #paths self.simplified_scattering_path = reader.read_param(params, 'simplified_scattering_path', '/host/input_data/simplified_scattering_params.txt') - #self.detailed_scattering_path = reader.read_param(params, 'detailed_scattering_path', None) + self.detailed_scatter_spectra_path = reader.read_param(params, 'path_to_detailed_scatter_spectra_dir', '.') self.efficiency_path = reader.read_param(params, 'efficiency_path', '') + self.final_states_file = reader.read_param(params, 'final_states_file', '') #options self.use_lineshape = reader.read_param(params, 'use_lineshape', True) self.detailed_or_simplified_lineshape = reader.read_param(params, 'detailed_or_simplified_lineshape', 'detailed') self.apply_efficiency = reader.read_param(params, 'apply_efficiency', False) self.return_frequency = reader.read_param(params, 'return_frequency', True) + self.molecular_final_states = reader.read_param(params, 'molecular_final_states', False) + # will be replaced with complex lineshape object if detailed lineshape is used + self.complexLineShape = None # get file content if needed # get efficiency dictionary @@ -125,16 +133,60 @@ def InternalConfigure(self, params): self.survival_prob, self.NScatters) elif self.lineshape=='detailed': - if not os.path.exists('./scatter_spectra_files'): - raise IOError('./scatter_spectra_files does not exist') + # check path exists + if 'scatter_spectra_file' in self.detailed_scatter_spectra_path: + full_path = self.detailed_scatter_spectra_path + self.detailed_scatter_spectra_path, _ = os.path.split(full_path) + else: + full_path = os.path.join(self.detailed_scatter_spectra_path, 'scatter_spectra_file') + + logger.info('Path to scatter_spectra_file: {}'.format(self.detailed_scatter_spectra_path)) + + + # lineshape params self.SimpParams = [self.scattering_sigma*2*math.sqrt(2*math.log(2)), self.survival_prob] + + # Setup and configure lineshape processor + complexLineShape_config = { + 'gases': ["H2","He"], + 'max_scatters': self.NScatters, + 'fix_scatter_proportion': True, + # When fix_scatter_proportion is True, set the scatter proportion for gas1 below + 'gas1_scatter_proportion': self.scatter_proportion, + # This is an important parameter which determines how finely resolved + # the scatter calculations are. 10000 seems to produce a stable fit with minimal slowdown, for ~4000 fake events. The parameter may need to + # be increased for larger datasets. + 'num_points_in_std_array': 10000, + 'B_field': self.B_field, + 'base_shape': 'dirac', + 'path_to_osc_strengths_files': self.detailed_scatter_spectra_path + } + logger.info('Setting up complex lineshape object') + self.complexLineShape = KrComplexLineShape("complexLineShape") + logger.info('Configuring complex lineshape') + self.complexLineShape.Configure(complexLineShape_config) + logger.info('Checking existence of scatter spectra files') + self.complexLineShape.check_existence_of_scatter_file() else: - raise ValueError("'detailed_or_simplified' is neither 'detailed' nor 'simplified'") + logger.error("'detailed_or_simplified' is neither 'detailed' nor 'simplified'") + return False else: self.lineshape = 'gaussian' self.SimpParams = [self.broadening] logger.info('Lineshape is Gaussian') + + # check final states file existence + if self.molecular_final_states: + logger.info('Going to use molecular final states from Bodine et al 2015') + with open(self.final_states_file, 'r') as infile: + a = json.load(infile) + index = np.where(np.array(a['Probability'])[:-1]>0) + self.final_state_array = [np.array(a['Binding energy'])[index], np.array(a['Probability'])[index]] + else: + logger.info('Not using molecular final state spectrum') + self.final_state_array = np.array([[0.], [1.]]) + return True @@ -142,7 +194,10 @@ def InternalConfigure(self, params): def InternalRun(self): if self.return_frequency: - ROIbound = [self.minf] + if self.maxf == None: + ROIbound = [self.minf] + else: + ROIbound = [self.minf, self.maxf] else: ROIbound = [self.Kmin, self.Kmax] @@ -213,10 +268,13 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4 if self.return_frequency: minf = ROIbound[0] - if efficiency_dict is not None: - maxf = max(efficiency_dict['frequencies']) + if len(ROIbound)==2: + maxf = ROIbound[1] else: - maxf = max(self.load_efficiency_curve()['frequencies']) + if efficiency_dict is not None: + maxf = max(efficiency_dict['frequencies']) + else: + maxf = max(self.load_efficiency_curve()['frequencies']) Kmax, Kmin = Energy(minf, B_field), Energy(maxf, B_field) else: Kmin, Kmax = ROIbound[0], ROIbound[1] @@ -226,10 +284,10 @@ def generate_unbinned_data(self, Q_mean, mass, ROIbound, S, B_1kev, nsteps=10**4 FWHM_convert = 2*math.sqrt(2*math.log(2)) if lineshape=='gaussian': max_energy = nstdevs*params[0] - min_energy = -1000 + min_energy = self.min_energy elif lineshape=='simplified_scattering' or lineshape=='simplified' or lineshape=='detailed_scattering' or lineshape=='detailed': max_energy = nstdevs/FWHM_convert*params[0] - min_energy = -1000 + min_energy = self.min_energy Kmax_eff = Kmax+max_energy #Maximum energy for data is slightly above Kmax>Q-m Kmin_eff = Kmin+min_energy #Minimum is slightly below Kmin 0.1: + logger.warning('Lineshape will shift spectrum by > 0.1 eV') + if min_x > 1.: + logger.warning('Lineshape will shift spectrum by > 1 eV') + raise ValueError('problem with std_eV_array()') return ans + # A gaussian centered at 0 eV with variable width, on the SELA def std_gaussian(self, sigma): x_array = self.std_eV_array() @@ -127,8 +148,7 @@ def normalize(self, f): def single_scatter_f(self, gas_type): energy_loss_array = self.std_eV_array() f = 0 * energy_loss_array - - input_filename = self.path_to_osc_strengths_files + gas_type + "OscillatorStrength.txt" + input_filename = os.path.join(self.path_to_osc_strengths_files, "{}OscillatorStrength.txt".format(gas_type)) energy_fOsc = ComplexLineShapeUtilities.read_oscillator_str_file(input_filename) fData = interpolate.interp1d(energy_fOsc[0], energy_fOsc[1], kind='linear') for i in range(len(energy_loss_array)): @@ -151,7 +171,7 @@ def another_scatter(self, input_spectrum, gas_type): return f_normed # Convolves the scatter functions and saves - # the results to a .npy file. + # the results to a .npy file. def generate_scatter_convolution_file(self): t = time.time() scatter_spectra_single_gas = {} @@ -179,7 +199,7 @@ def generate_scatter_convolution_file(self): total_scatter = self.normalize(signal.convolve(H2_scatter, Kr_scatter, mode='same')) scatter_spectra['{}_{}'.format(self.gases[0], self.gases[1])]['{}_{}'.format(str(j).zfill(2), str(i-j).zfill(2))] = total_scatter np.save( - self.path_to_osc_strengths_files+'scatter_spectra_file/scatter_spectra.npy', + os.path.join(self.path_to_osc_strengths_files, 'scatter_spectra_file/scatter_spectra.npy'), scatter_spectra ) elapsed = time.time() - t @@ -191,17 +211,17 @@ def generate_scatter_convolution_file(self): # If not, this function calls generate_scatter_convolution_file. # This function also checks to make sure that the scatter file have the correct # number of entries and correct number of points in the SELA, and if not, it generates a fresh file. - # When the variable regenerate is set as True, it generates a fresh file + # When the variable regenerate is set as True, it generates a fresh file def check_existence_of_scatter_file(self, regenerate = False): gases = self.gases if regenerate == True: logger.info('generate fresh scatter file') self.generate_scatter_convolution_file() - else: + else: stuff_in_dir = os.listdir(self.path_to_osc_strengths_files) if 'scatter_spectra_file' not in stuff_in_dir: logger.info('Scatter spectra folder not found, generating') - os.mkdir(self.path_to_osc_strengths_files+'scatter_spectra_file') + os.mkdir(os.path.join(self.path_to_osc_strengths_files,'scatter_spectra_file')) time.sleep(2) self.generate_scatter_convolution_file() else: @@ -209,7 +229,7 @@ def check_existence_of_scatter_file(self, regenerate = False): strippeddirs = [s.strip('\n') for s in directory] if 'scatter_spectra.npy' not in strippeddirs: self.generate_scatter_convolution_file() - test_file = self.path_to_osc_strengths_files+'scatter_spectra_file/scatter_spectra.npy' + test_file = os.path.join(self.path_to_osc_strengths_files, 'scatter_spectra_file/scatter_spectra.npy') test_dict = np.load(test_file, allow_pickle = True) if list(test_dict.item().keys())[0] != '{}_{}'.format(gases[0], gases[1]): logger.info('first entry not matching, generating fresh files') @@ -249,6 +269,8 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt current_working_spectrum = self.std_lorenztian_17keV() elif emitted_peak == 'shake': current_working_spectrum = self.shakeSpectrumClassInstance.shake_spectrum() + elif emitted_peak == 'dirac': + current_working_spectrum = self.std_dirac() current_working_spectrum = self.convolve_gaussian(current_working_spectrum, gauss_FWHM_eV) zeroth_order_peak = current_working_spectrum current_full_spectrum += current_working_spectrum @@ -276,6 +298,7 @@ def make_spectrum(self, gauss_FWHM_eV, prob_parameter, scatter_proportion, emitt current_working_spectrum = \ self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same')) current_full_spectrum += current_working_spectrum*(prob_parameter*q)**(n) + return current_full_spectrum # Produces a spectrum in real energy that can now be evaluated off of the SELA. @@ -294,12 +317,12 @@ def spectrum_func(self, x_keV, *p0): amplitude = p0[2] prob_parameter = p0[3] scatter_proportion = p0[4] - + line_pos_eV = line_pos_keV*1000. x_eV_minus_line = x_eV - line_pos_eV zero_idx = np.r_[np.where(x_eV_minus_line< en_loss_array_min)[0],np.where(x_eV_minus_line>en_loss_array_max)[0]] nonzero_idx = [i for i in range(len(x_keV)) if i not in zero_idx] - + full_spectrum = self.make_spectrum(FWHM_G_eV, prob_parameter, scatter_proportion) full_spectrum_rev = ComplexLineShapeUtilities.flip_array(full_spectrum) f_intermediate[nonzero_idx] = np.interp(x_eV_minus_line[nonzero_idx],en_array_rev,full_spectrum_rev) @@ -340,8 +363,8 @@ def fit_data(self, freq_bins, data_hist_freq, print_params=True): amplitude_guess = np.sum(data_hist)/2 prob_parameter_guess = 0.5 scatter_proportion_guess = 0.9 - p0_guess = [FWHM_guess, line_pos_guess, amplitude_guess, prob_parameter_guess, scatter_proportion_guess] - p0_bounds = ([FWHM_eV_min, line_pos_keV_min, amplitude_min, prob_parameter_min, scatter_proportion_min], + p0_guess = [FWHM_guess, line_pos_guess, amplitude_guess, prob_parameter_guess, scatter_proportion_guess] + p0_bounds = ([FWHM_eV_min, line_pos_keV_min, amplitude_min, prob_parameter_min, scatter_proportion_min], [FWHM_eV_max, line_pos_keV_max, amplitude_max, prob_parameter_max, scatter_proportion_max]) # Actually do the fitting params , cov = curve_fit(self.spectrum_func, bins_keV_nonzero, data_hist_nonzero, sigma=data_hist_err, p0=p0_guess, bounds=p0_bounds) @@ -362,7 +385,7 @@ def fit_data(self, freq_bins, data_hist_freq, print_params=True): prob_parameter_fit_err = perr[3] scatter_proportion_fit_err = perr[4] total_counts_fit_err = amplitude_fit_err - + fit = self.spectrum_func(bins_keV,*params) line_pos_Hz_fit , line_pos_Hz_fit_err = ComplexLineShapeUtilities.energy_guess_to_frequency(line_pos_keV_fit, line_pos_keV_fit_err, self.B_field) @@ -421,7 +444,7 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'): p = self.scatter_proportion q = 1 - p scatter_spectra = np.load( - current_path + 'scatter_spectra_file/scatter_spectra.npy', allow_pickle = True + os.path.join(current_path,'scatter_spectra_file/scatter_spectra.npy'), allow_pickle = True ) en_array = self.std_eV_array() current_full_spectrum = np.zeros(len(en_array)) @@ -429,6 +452,8 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'): current_working_spectrum = self.std_lorenztian_17keV() elif emitted_peak == 'shake': current_working_spectrum = self.shakeSpectrumClassInstance.shake_spectrum() + elif emitted_peak == 'dirac': + current_working_spectrum = self.std_dirac() current_working_spectrum = self.convolve_gaussian(current_working_spectrum, gauss_FWHM_eV) zeroth_order_peak = current_working_spectrum current_full_spectrum += current_working_spectrum @@ -455,6 +480,7 @@ def make_spectrum_1(self, gauss_FWHM_eV, prob_parameter, emitted_peak='shake'): current_working_spectrum = \ self.normalize(signal.convolve(zeroth_order_peak, current_working_spectrum, mode='same')) current_full_spectrum += current_working_spectrum*(prob_parameter*q)**(n) + return current_full_spectrum def spectrum_func_1(self, x_keV, *p0): @@ -476,7 +502,7 @@ def spectrum_func_1(self, x_keV, *p0): zero_idx = np.r_[np.where(x_eV_minus_line< en_loss_array_min)[0],np.where(x_eV_minus_line>en_loss_array_max)[0]] nonzero_idx = [i for i in range(len(x_keV)) if i not in zero_idx] - full_spectrum = self.make_spectrum_1(FWHM_G_eV, prob_parameter,) + full_spectrum = self.make_spectrum_1(FWHM_G_eV, prob_parameter,emitted_peak=self.base_shape) full_spectrum_rev = ComplexLineShapeUtilities.flip_array(full_spectrum) f_intermediate[nonzero_idx] = np.interp(x_eV_minus_line[nonzero_idx],en_array_rev,full_spectrum_rev) f[nonzero_idx] += amplitude*f_intermediate[nonzero_idx]/np.sum(f_intermediate[nonzero_idx]) @@ -506,8 +532,8 @@ def fit_data_1(self, freq_bins, data_hist_freq): line_pos_guess = bins_keV[np.argmax(data_hist)] amplitude_guess = np.sum(data_hist)/2 prob_parameter_guess = 0.5 - p0_guess = [FWHM_guess, line_pos_guess, amplitude_guess, prob_parameter_guess] - p0_bounds = ([FWHM_eV_min, line_pos_keV_min, amplitude_min, prob_parameter_min], + p0_guess = [FWHM_guess, line_pos_guess, amplitude_guess, prob_parameter_guess] + p0_bounds = ([FWHM_eV_min, line_pos_keV_min, amplitude_min, prob_parameter_min], [FWHM_eV_max, line_pos_keV_max, amplitude_max, prob_parameter_max]) # Actually do the fitting params , cov = curve_fit(self.spectrum_func_1, bins_keV_nonzero, data_hist_nonzero, sigma=data_hist_err, p0=p0_guess, bounds=p0_bounds) @@ -526,7 +552,7 @@ def fit_data_1(self, freq_bins, data_hist_freq): amplitude_fit_err = perr[2] prob_parameter_fit_err = perr[3] total_counts_fit_err = amplitude_fit_err - + fit = self.spectrum_func_1(bins_keV,*params) line_pos_Hz_fit , line_pos_Hz_fit_err = ComplexLineShapeUtilities.energy_guess_to_frequency(line_pos_keV_fit, line_pos_keV_fit_err, self.B_field) @@ -568,4 +594,4 @@ def fit_data_1(self, freq_bins, data_hist_freq): 'amplitude_fit_err': amplitude_fit_err, 'data_hist_freq': data_hist_freq } - return dictionary_of_fit_results \ No newline at end of file + return dictionary_of_fit_results diff --git a/setup.py b/setup.py index 56aa86ef..170d5474 100644 --- a/setup.py +++ b/setup.py @@ -20,7 +20,7 @@ requirements = [] extras_require = { - 'core':['colorlog'], + 'core':['colorlog', 'iminuit'], 'doc': ['sphinx','sphinx_rtd_theme','sphinxcontrib-programoutput', 'six', 'colorlog'] } diff --git a/tests/Fake_data_generator_test.py b/tests/Fake_data_generator_test.py index 22a41025..6e729297 100644 --- a/tests/Fake_data_generator_test.py +++ b/tests/Fake_data_generator_test.py @@ -20,16 +20,20 @@ def test_data_generation(self): "apply_efficiency": False, "efficiency_path": "./combined_energy_corrected_eff_at_quad_trap_frequencies.json", "simplified_lineshape_path": None, + "path_to_detailed_scatter_spectra_dir": '/host', "detailed_or_simplified_lineshape": "detailed", #"simplified" or "detailed" "use_lineshape": False, # if False only gaussian smearing is applied "return_frequency": True, "scattering_sigma": 18.6, # only used if use_lineshape = True - "scattering_prob": 0.77, # only used if use_lineshape = True + "survival_prob": 0.77, # only used if use_lineshape = True + "scatter_proportion": 0.8, # only used if use_lineshape = True and lineshape = detailed "B_field": 0.9578186017836624, "S": 4500, # number of tritium events - "n_steps": 1000, # stepsize for pseudo continuous data is: (Kmax_eff-Kmin_eff)/nsteps + "n_steps": 10000, # stepsize for pseudo continuous data is: (Kmax_eff-Kmin_eff)/nsteps "A_b": 1e-10, # background rate 1/eV/s - "poisson_stats": True + "poisson_stats": True, + "molecular_final_states": False, + "final_states_file": "../mermithid/misc/saenz_mfs.json" } specGen = FakeDataGenerator("specGen") diff --git a/tests/Fitters_test.py b/tests/Fitters_test.py new file mode 100644 index 00000000..7029d18b --- /dev/null +++ b/tests/Fitters_test.py @@ -0,0 +1,142 @@ +""" +Script to test fit processors +Author: C. Claessens +Date: August 16, 2020 + +""" + +import unittest + +from morpho.utilities import morphologging, parser +logger = morphologging.getLogger(__name__) + +import matplotlib.pyplot as plt +import numpy as np + +class FittersTest(unittest.TestCase): + + def test_iminuit(self): + logger.info('iMinuit test') + from iminuit import Minuit + + def f(x, y, z): + return (x - 2) ** 2 + (y - 3) ** 2 + (z - 4) ** 2 + + def g(params): + return f(*params) + + m = Minuit(f, x=0, y=0, z=0) + + m.migrad() # run optimiser + print(m.values) # {'x': 2,'y': 3,'z': 4} + + m.hesse() # run covariance estimator + print(m.errors) # {'x': 1,'y': 1,'z': 1} + + # repeat using from_array_func + m2 = Minuit(g, [0, 0, 0], name=['a1', 'b1', 'c1']) + m2.errors = [0.1, 0.1, 0.1] + m2.errordef = 1 + m2.print_level = 0 + m2.migrad() + print(m2.values) + + logger.info('iMinuit test done') + + def test_BinnedDataFitter(self): + from mermithid.processors.Fitters import BinnedDataFitter + + logger.info('iMinuit fit test') + config_dict = { + 'variables': 'N', + 'bins': np.linspace(-3, 3, 100), + 'parameter_names': ['A', 'mu', 'sigma'], + 'initial_values': [100, 0, 1], + 'limits': [[0, None], [None, None], [0, None]], + 'constrained_parameter_indices': [], + 'constrained_parameter_means': [0.5], + 'constrained_parameter_widths': [1], + 'binned_data': True, + 'print_level': 0 + + } + + random_data = {'K': np.random.randn(10000)*0.5+1} + + # histogramming could be done by processor + binned_data, _ = np.histogram(random_data['K'], config_dict['bins']) + binned_data_dict = {'N': binned_data} + + # setup processor + negll_fitter = BinnedDataFitter('iminuit_processor') + negll_fitter.Configure(config_dict) + + # define new model and overwrite processor's model + def gaussian(x, A, mu, sigma): + """ + This is the same function that is implemented as default model + """ + f = A*(1/(sigma*np.sqrt(2*np.pi)))*np.exp(-(((x-mu)/sigma)**2.)/2.) + return f + + negll_fitter.model = gaussian + + # hand over data and run + negll_fitter.data = binned_data_dict #random_data + negll_fitter.Run() + + # collect fit results + results = negll_fitter.results + result_list = results['param_values'] + error_list = results['param_errors'] + + for k in results.keys(): + logger.info('{}: {}'.format(k, results[k])) + + + + # get bins histogram and fit curve from processor for plotting + x = negll_fitter.bin_centers + hist = negll_fitter.hist + hist_fit = negll_fitter.model(x, *result_list) + + # calculate normalized residuals + residuals = (hist-hist_fit)/np.sqrt(hist_fit) + + plt.rcParams.update({'font.size': 20}) + plt.figure(figsize=(10,10)) + plt.subplot(211) + plt.errorbar(x, hist, np.sqrt(hist), drawstyle='steps-mid', label='Binned data') + plt.plot(x, negll_fitter.model(x, *result_list), label='Fit') + plt.xlabel(negll_fitter.namedata) + plt.ylabel('Counts') + plt.legend() + + plt.subplot(212) + plt.scatter(x, residuals, label='Pearson residuals') + plt.axhline(np.nanmean(residuals)) + plt.xlabel(negll_fitter.namedata) + plt.ylabel('Residuals') + plt.legend() + + + plt.savefig('iminuit_fit.png') + + + + +if __name__ == '__main__': + + args = parser.parse_args(False) + + + logger = morphologging.getLogger('morpho', + level=args.verbosity, + stderr_lb=args.stderr_verbosity, + propagate=False) + logger = morphologging.getLogger(__name__, + level=args.verbosity, + stderr_lb=args.stderr_verbosity, + propagate=False) + + unittest.main() \ No newline at end of file diff --git a/tests/IO_test.py b/tests/IO_test.py index 7a1050b4..0ea2cd39 100644 --- a/tests/IO_test.py +++ b/tests/IO_test.py @@ -75,9 +75,9 @@ def test_MutliChannelReader(self): plt.ylabel('N') plt.subplot(212) - plt.hist(data['a']['TrueStartFrequenciesCut']-24.5e9, bins=bins-24.5e9, label='channel a: {} counts'.format(len(data['a']['TrueStartFrequenciesCut']))) - plt.hist(data['b']['TrueStartFrequenciesCut']-24.5e9, bins=bins-24.5e9, label='channel b: {} counts'.format(len(data['b']['TrueStartFrequenciesCut']))) - plt.hist(data['c']['TrueStartFrequenciesCut']-24.5e9, bins=bins-24.5e9, label='channel c: {} counts'.format(len(data['c']['TrueStartFrequenciesCut']))) + plt.hist(np.array(data['a']['TrueStartFrequenciesCut'])-24.5e9, bins=bins-24.5e9, label='channel a: {} counts'.format(len(data['a']['TrueStartFrequenciesCut']))) + plt.hist(np.array(data['b']['TrueStartFrequenciesCut'])-24.5e9, bins=bins-24.5e9, label='channel b: {} counts'.format(len(data['b']['TrueStartFrequenciesCut']))) + plt.hist(np.array(data['c']['TrueStartFrequenciesCut'])-24.5e9, bins=bins-24.5e9, label='channel c: {} counts'.format(len(data['c']['TrueStartFrequenciesCut']))) plt.xlabel('Start frequencies') plt.ylabel('N') plt.legend() diff --git a/tests/Misc_test.py b/tests/Misc_test.py index 54e0edcb..d1f05019 100644 --- a/tests/Misc_test.py +++ b/tests/Misc_test.py @@ -8,27 +8,50 @@ import unittest -from morpho.utilities import morphologging, read_param +from morpho.utilities import morphologging, parser logger = morphologging.getLogger(__name__) +import matplotlib.pyplot as plt +import numpy as np -class TritiumTests(unittest.TestCase): - from mermithid.processors.misc import FrequencyEnergyConversionProcessor - - freq_data = [27.9925*10**9, 27.0094*10**9, - 26.4195*10**9, 26.4169*10**9, - 26.3460*10**9, 26.3457*10**9] - logger.info("Will convert the following frequencies: %s"%freq_data) - logger.debug("At 1 T, these correspond to kinetic energies (in keV) of " + - "[0, 18.6, 30.424, 30.477, 31.934, 31.942]") - - freq_energy_dict = { - "B": 1 - } - - freq_proc = FrequencyEnergyConversionProcessor("freq_energy_processor") - freq_proc.Configure(freq_energy_dict) - freq_proc.frequencies = freq_data - freq_proc.Run() - - logger.info("Resulting energies: %s"%freq_proc.energies) +class MiscTest(unittest.TestCase): + + def test_FreqConversionTest(self): + from mermithid.processors.misc import FrequencyEnergyConversionProcessor + + freq_data = [27.9925*10**9, 27.0094*10**9, + 26.4195*10**9, 26.4169*10**9, + 26.3460*10**9, 26.3457*10**9] + logger.info("Will convert the following frequencies: %s"%freq_data) + logger.debug("At 1 T, these correspond to kinetic energies (in keV) of " + + "[0, 18.6, 30.424, 30.477, 31.934, 31.942]") + + freq_energy_dict = { + "B": 1 + } + + freq_proc = FrequencyEnergyConversionProcessor("freq_energy_processor") + freq_proc.Configure(freq_energy_dict) + freq_proc.frequencies = freq_data + freq_proc.Run() + + logger.info("Resulting energies: %s"%freq_proc.energies) + + + + +if __name__ == '__main__': + + args = parser.parse_args(False) + + + logger = morphologging.getLogger('morpho', + level=args.verbosity, + stderr_lb=args.stderr_verbosity, + propagate=False) + logger = morphologging.getLogger(__name__, + level=args.verbosity, + stderr_lb=args.stderr_verbosity, + propagate=False) + + unittest.main() \ No newline at end of file diff --git a/tests/test.sh b/tests/test.sh index 85c649e3..8483f480 100644 --- a/tests/test.sh +++ b/tests/test.sh @@ -5,3 +5,4 @@ python Misc_test.py python Tritium_test.py python TritiumAndEfficiencyBinner_test.py python Fake_data_generator_test.py +python Fitters_test.py