diff --git a/_test/test_gen_sqw_workflow/test_gen_sqw_accumulate_sqw_tests_nxspe.m b/_test/test_gen_sqw_workflow/test_gen_sqw_accumulate_sqw_tests_nxspe.m new file mode 100644 index 0000000000..c7d8dafd81 --- /dev/null +++ b/_test/test_gen_sqw_workflow/test_gen_sqw_accumulate_sqw_tests_nxspe.m @@ -0,0 +1,863 @@ +classdef test_gen_sqw_accumulate_sqw_tests_nxspe < TestCaseWithSave + % Series of tests of gen_sqw and associated functions + % generated using multiple Matlab workers. + % + % Optionally writes results to output file to compare with previously + % saved sample test results + %--------------------------------------------------------------------- + % Usage: + % + %1) Normal usage: + % Run all unit tests and compare their results with previously saved + % results stored in test_gen_sqw_accumulate_sqw_output.mat file + % located in the same folder as this function: + % + %>>runtests test_gen_sqw_accumulate_sqw_sep_session + %--------------------------------------------------------------------- + %2) Run particular test case from the suite: + % + %>>tc = test_gen_sqw_accumulate_sqw_sep_session(); + %>>tc.test_[particular_test_name] e.g.: + %>>tc.test_accumulate_sqw14(); + %or + %>>tc.test_gen_sqw(); + %--------------------------------------------------------------------- + %3) Generate test file to store test results to compare with them later + % (it stores test results into tmp folder.) + % + %>>tc=test_gen_sqw_accumulate_sqw_sep_session('save'); + %>>tc.save(): + properties + % properties to use as input for data + test_data_path; + test_functions_path; + par_file; + nfiles_max=6; + + pars; + scale; + + proj; + gen_sqw_par={}; + % files; + spe_file={[]}; + + instrum + sample + % + % the field stores initial configuration, was in place when test + % was started to run + initial_config; + % + % the property describes common name of the test files allowing to + % distinguish these files from the files, generated by other type + % of test + test_pref = 'nomex'; + + working_dir + + nxspedir + nxspefile + + end + + methods(Static) + function new_names = rename_file_list(input_list,new_ext) + % change extension for list of files + if ~iscell(input_list) + input_list = {input_list}; + end + new_names = cell(1,numel(input_list)); + for i=1:numel(input_list) + fls = input_list{i}; + [fpath,fn,~] = fileparts(fls); + flt = fullfile(fpath,[fn,new_ext]); + new_names{i} = flt; + if is_file(fls) + movefile(fls,flt,'f'); + end + end + end + end + + methods + function obj=test_gen_sqw_accumulate_sqw_tests_nxspe(name) + obj.nxspedir = 'C:\Users\nvl96446\STFC\PACE\nxspe\merlin'; %maps'; + obj.nxspefile = [obj.nxspedir '\MER62983_13.2meV_1to1.nxspe']; %MAP45862_250meV_1to1.nxspe']; + if ~exist(obj.nxspefile, 'file') + disp([obj.nxspefile ' not there']); + else + disp('it''s there'); + end + + obj.nfiles_max=1; + en=cell(1,obj.nfiles_max); + efix=zeros(1,obj.nfiles_max); + psi=zeros(1,obj.nfiles_max); + omega=zeros(1,obj.nfiles_max); + dpsi=zeros(1,obj.nfiles_max); + gl=zeros(1,obj.nfiles_max); + gs=zeros(1,obj.nfiles_max); + for i=1:obj.nfiles_max + efix(i)=230+0.5*i; % different ei for each file + en{i}=0.05*efix(i):0.2+i/50:0.95*efix(i); % different energy bins for each file + psi(i)=90-i+1; + omega(i)=10+i/2; + dpsi(i)=0.1+i/10; + gl(i)=3-i/6; + gs(i)=2.4+i/7; + end + psi=90:-1:90-obj.nfiles_max+1; + + emode=1; + alatt=[4.4,5.5,6.6]; + angdeg=[100,105,110]; + u=[1.02,0.99,0.02]; + v=[0.025,-0.01,1.04]; + + obj.gen_sqw_par={en,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs}; + end + + function test_dummy(obj,varargin) + end + + function test_gen_nxspe(obj,varargin) + extra_path = '..\test_TF_components'; + addpath(extra_path); + select = 1:1; + en =obj.gen_sqw_par{1}(select); + efix=obj.gen_sqw_par{2}(select); + emode=obj.gen_sqw_par{3}; + alatt=obj.gen_sqw_par{4}; + angdeg=obj.gen_sqw_par{5}; + u=obj.gen_sqw_par{6}; + v=obj.gen_sqw_par{7}; + psi=obj.gen_sqw_par{8}(select); + omega=obj.gen_sqw_par{9}(select); + dpsi=obj.gen_sqw_par{10}(select); + gl=obj.gen_sqw_par{11}(select); + gs=obj.gen_sqw_par{12}(select); + + [tmp_files,grid2,pix_data_range2]=gen_sqw (obj.nxspefile, ... + '', 'outfile', efix, emode, alatt, angdeg,... + u, v, psi, omega, ... + dpsi, gl, gs); + end + + %{ + function obj=gen_sqw_accumulate_sqw_tests_common(test_class_name,test_prefix) + % The constructor for class, which is the common part of all + % MPI-based gen_sqw system tests. + % + % Should be used as + % + % >> test_gen_sqw_accumulate_sqw % Compares with + % previously saved results in + % test_gen_sqw_accumulate_sqw_output.mat + % % in the same + % folder as this + % function + % >> test_gen_sqw_accumulate_sqw ('save') % Save to + % test_gen_sqw_accumulate_sqw_output.mat + % + % Reads previously created test data sets. + + + obj = obj@TestCaseWithSave(test_class_name,fullfile(fileparts(mfilename('fullpath')),'test_gen_sqw_accumulate_sqw_output.mat')); + if exist('test_prefix', 'var') + obj.test_pref = test_prefix; + end + % do other initialization + obj.comparison_par={ 'min_denominator', 0.01, 'ignore_str', 1}; + obj.tol = 1.e-5; + + pths = horace_paths; + obj.test_functions_path = pths.test_common_func; + + addpath(obj.test_functions_path); + + if isempty(obj.new_working_folder) + hc = hor_config; + obj.working_dir = hc.working_directory; + else + obj.working_dir = obj.new_working_folder; + end + + if ~is_folder(obj.working_dir) + mkdir(obj.working_dir); + end + + % build test file names + obj.spe_file=cell(1,obj.nfiles_max); + wkd = obj.working_dir; + for i=1:obj.nfiles_max + obj.spe_file{i}=fullfile(wkd ,... + ['gen_sqw_acc_sqw_spe_',test_prefix,num2str(i),'.nxspe']); + end + + + common_data = fullfile(fileparts(fileparts(mfilename('fullpath'))),'common_data'); + %this.par_file=fullfile(this.results_path,'96dets.par'); + obj.par_file=fullfile(common_data,'gen_sqw_96dets.nxspe'); + + + % initiate test parameters + en=cell(1,obj.nfiles_max); + efix=zeros(1,obj.nfiles_max); + psi=zeros(1,obj.nfiles_max); + omega=zeros(1,obj.nfiles_max); + dpsi=zeros(1,obj.nfiles_max); + gl=zeros(1,obj.nfiles_max); + gs=zeros(1,obj.nfiles_max); + for i=1:obj.nfiles_max + efix(i)=35+0.5*i; % different ei for each file + en{i}=0.05*efix(i):0.2+i/50:0.95*efix(i); % different energy bins for each file + psi(i)=90-i+1; + omega(i)=10+i/2; + dpsi(i)=0.1+i/10; + gl(i)=3-i/6; + gs(i)=2.4+i/7; + end + psi=90:-1:90-obj.nfiles_max+1; + + emode=1; + alatt=[4.4,5.5,6.6]; + angdeg=[100,105,110]; + u=[1.02,0.99,0.02]; + v=[0.025,-0.01,1.04]; + + obj.gen_sqw_par={en,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs}; + + obj.pars=[1000,8,2,4,0]; % [Seff,SJ,gap,gamma,bkconst] + obj.scale=0.3; + % build test files if they have not been build + obj=build_test_files(obj); + + obj.save(); + end + %} + % + % + function [en,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs]=unpack(this,varargin) + if nargin>1 + n_elem = varargin{1}; + if numel(n_elem) >1 + select = n_elem; + else + select = 1:n_elem; + end + else + n_elem = numel(this.gen_sqw_par{1}); + select = 1:n_elem; + end + en =this.gen_sqw_par{1}(select); + efix=this.gen_sqw_par{2}(select); + emode=this.gen_sqw_par{3}; + alatt=this.gen_sqw_par{4}; + angdeg=this.gen_sqw_par{5}; + u=this.gen_sqw_par{6}; + v=this.gen_sqw_par{7}; + psi=this.gen_sqw_par{8}(select); + omega=this.gen_sqw_par{9}(select); + dpsi=this.gen_sqw_par{10}(select); + gl=this.gen_sqw_par{11}(select); + gs=this.gen_sqw_par{12}(select); + end + % + function obj=build_test_files(obj,spe_files) + if ~exist('spe_files','var') + spe_files = obj.spe_file; + end + + % ===================================================================================================================== + % Make instrument and sample + % ===================================================================================================================== + wmod=IX_moderator('AP2',12,35,'ikcarp',[3,25,0.3],'',[],0.12,0.12,0.05,300); + wap=IX_aperture(-2,0.067,0.067); + wchop=IX_fermi_chopper(1.8,600,0.1,1.3,0.003); + instrument_ref.moderator=wmod; + instrument_ref.aperture=wap; + instrument_ref.fermi_chopper=wchop; + sample_ref=IX_sample('PCSMO',true,[1,1,0],[0,0,1],'cuboid',[0.04,0.05,0.02],1.6,300); + + instrument=repmat(instrument_ref,1,numel(spe_files)); + for i=1:numel(instrument) + instrument(i).IX_fermi_chopper.frequency=100*i; + end + obj.instrum = IX_inst_DGfermi(instrument_ref); + obj.sample = sample_ref; + + + + file_exist = cellfun(@(fn)(exist(fn,'file') == 2),spe_files); + if all(file_exist) + obj.add_to_files_cleanList(spe_files{:}); + return; + end + spe_files = spe_files(~file_exist); + n_files = numel(spe_files); + + + %sample_1=sample_ref; + %sample_2=sample_ref; + %sample_2.temperature=350; + + + % ===================================================================================================================== + % Make spe files + % ===================================================================================================================== + % for the purposes of consistency, test files have to be always + % generated by single thread, as multithreading causes pixels + % permutation within a bin, and then random function adds + % various amounts of noise to various detectors according to the + % ordering + + % uses Matlab single thread to generate source sqw files + clob_hc = set_temporary_config_options(hor_config, 'use_mex', false); + clob_pc = set_temporary_config_options(parallel_config, 'threads', 1); + + [en,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs]=unpack(obj); + for i=1:n_files + if ~is_file(spe_files{i}) + simulate_spe_testfunc (en{i}, obj.par_file,spe_files{i}, @sqw_sc_hfm_testfunc, obj.pars, obj.scale,... + efix(i), emode, alatt, angdeg, u, v, psi(i), omega(i), dpsi(i), gl(i), gs(i)); + end + end + + obj.add_to_files_cleanList(spe_files{:}); + end + % + function test_gen_sqw(obj,varargin) + %------------------------------------------------------------- + if obj.skip_test + skipTest(fprintf('test_gen_sqw_%s is disabled',obj.test_pref)); + end + if nargin> 1 + % running in single test method mode. + obj.setUp(); + co1 = onCleanup(@()obj.tearDown()); + end + %------------------------------------------------------------- + + + % build test files if they have not been build + obj=build_test_files(obj); + % generate the names of the output sqw files + + sqw_file=cell(1,obj.nfiles_max); + file_pref = obj.test_pref; + wkdir = obj.working_dir; + for i=1:obj.nfiles_max + sqw_file{i}=fullfile(wkdir ,['test_gen_sqw_',file_pref ,num2str(i),'.sqw']); % output sqw file + end + + sqw_file_123456=fullfile(wkdir ,['sqw_123456_',file_pref,'.sqw']); % output sqw file + sqw_file_145623=fullfile(wkdir ,['sqw_145623_',file_pref,'.sqw']); % output sqw file + if ~obj.save_output + cleanup_obj1=onCleanup(@()obj.delete_files(sqw_file_123456,sqw_file_145623,sqw_file{:})); + end + % --------------------------------------- + % Test gen_sqw --------------------------------------- + clConf = set_temporary_config_options('hor_config','delete_tmp',true); + + [en,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs]=unpack(obj); + %hc.threads = 1; + + + [tmp_files,grid1,pix_data_range1]=gen_sqw (obj.spe_file, '', ... + sqw_file_123456, efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs); + for i=1:numel(tmp_files) + assertFalse(is_file(tmp_files{i})); + end + %hc.build_sqw_in_parallel=0; + [tmp_files,grid2,pix_data_range2]=gen_sqw (obj.spe_file([1,4,5,6,2,3]),... + '', sqw_file_145623, efix([1,4,5,6,2,3]), emode, alatt, angdeg,... + u, v, psi([1,4,5,6,2,3]), omega([1,4,5,6,2,3]), ... + dpsi([1,4,5,6,2,3]), gl([1,4,5,6,2,3]), gs([1,4,5,6,2,3])); + for i=1:numel(tmp_files) + assertFalse(is_file(tmp_files{i})); + end + + assertEqual(grid1,grid2); + assertElementsAlmostEqual(pix_data_range1,pix_data_range2); + + % Make some cuts: --------------- + obj.proj.u=[1,0,0.1]; obj.proj.v=[0,0,1]; + + + % Check cuts from gen_sqw output with spe files in a different + % order are the same + [ok,mess,dummy_w1,w1b]=is_cut_equal(sqw_file_123456,sqw_file_145623,obj.proj,[-1.5,0.025,0],[-2.1,-1.9],[-0.5,0.5],[-Inf,Inf]); + assertTrue(ok,['Cuts from gen_sqw output with spe files in a different order are not the same: ',mess]); + % Test against saved or store to save later + obj.assertEqualToTolWithSave(w1b,'ignore_str',true,'tol',1.e-7); + + + w1a=cut_sqw(sqw_file_123456,obj.proj,[-1.5,0.025,0],[-2.1,-1.9],[-0.5,0.5],[-Inf,Inf]); + % Test against saved or store to save later + obj.assertEqualToTolWithSave(w1a,'ignore_str',true,'tol',1.e-7); + end + % + function test_gen_sqw_sym(obj,varargin) + %------------------------------------------------------------- + if obj.save_output + return; + end + if obj.skip_test + skipTest(fprintf('test_gen_sqw_sym_%s is disabled',obj.test_pref)); + end + if nargin> 1 + % running in single test method mode. + obj.setUp(); + co1 = onCleanup(@()obj.tearDown()); + end + %------------------------------------------------------------- + + + % build test files if they have not been build + obj=build_test_files(obj); + % generate the names of the output sqw files + + file_pref = obj.test_pref; + wkdir = obj.working_dir; + + + sqw_file_base=fullfile(wkdir ,['sqw_sym_base_',file_pref,'.sqw']); % output basic sqw file + sqw_file_sym =fullfile(wkdir ,['sqw_sym_reflected_',file_pref,'.sqw']); % symmetrized sqw file + if ~obj.save_output + cleanup_obj1=onCleanup(@()obj.delete_files(sqw_file_base,sqw_file_sym)); + end + % --------------------------------------- + + [~,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs]=unpack(obj); + + hpc_config_cleanup = set_temporary_config_options(hpc_config, ... + 'build_sqw_in_parallel', false, ... + 'combine_sqw_using', 'mex_code' ... + ); + hc_config_cleanup = set_temporary_config_options(hor_config, 'delete_tmp', true); + + % Test symmetrisation --------------------------------------- + % Standard reference file. Done serially with mex combining for + % speed. + gen_sqw (obj.spe_file, '', sqw_file_base,... + efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs); %,... + + % symmetrise in memory + w_inm = read_sqw(sqw_file_base); + sym = SymopReflection([0,1,0], [0,0,1], [0,0,0]); + + w_mem_sym=symmetrise_sqw(w_inm, sym); + % return the configuration to the state, + % specified by tests + + % check all pixels are there nothing have been lost at + % symmetrisation + assertEqual(w_inm.pix.num_pixels,w_mem_sym.pix.num_pixels) + + gen_sqw(obj.spe_file, '', sqw_file_sym,... + efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs,... + 'transform_sqw',@(x)symmetrise_sqw(x,sym)); + + + loc_proj=struct('u',u,'v',v); + + w1_inf_sym = cut(sqw_file_sym,loc_proj,[-1.5,0.025,0],[-2.1,-1.9],[-0.5,0.5],[-Inf,Inf]); + w1_inm_sym = cut(w_mem_sym,loc_proj,[-1.5,0.025,0],[-2.1,-1.9],[-0.5,0.5],[-Inf,Inf]); + % Uncomment to see the cut shapes + % acolor('r') + % plot(w1_inf_sym) + % acolor('k') + % pd(w1_inm_sym) + + assertEqualToTol(w1_inf_sym,w1_inm_sym,'ignore_str',true,'tol',1.e-6) + + w_inf_sym = read_sqw(sqw_file_sym); + assertEqual(w_inf_sym.npixels,w_mem_sym.npixels) + %TODO: with proper projection the error should be lower. The + %difference due to various borders added in different places + clear hc_config_cleanup; + skipTest('Memory based and file-based symmetrisation ranges are different Ticket #798'); + assertEqualToTol(w_mem_sym,w_inf_sym,'ignore_str',true,'tol',2.e-5) + + + + end + + function test_accumulate_sqw14(obj,varargin) + %------------------------------------------------------------- + if obj.skip_test + skipTest(fprintf('test_accumulate_sqw14_%s is disabled',obj.test_pref)); + end + if nargin> 1 + % running in single test method mode. + obj.setUp(); + clobS = onCleanup(@()obj.tearDown()); + end + %------------------------------------------------------------- + clConf = set_temporary_config_options('hor_config','delete_tmp',true); + file_pref = obj.test_pref; + wk_dir = obj.working_dir; + + sqw_file_accum=fullfile(wk_dir ,['test_sqw_accum_sqw14_',file_pref,'.sqw']); + sqw_file_14=fullfile(wk_dir ,['test_sqw_14_',file_pref,'.sqw']); + clobR=onCleanup(@()obj.delete_files(sqw_file_14,sqw_file_accum)); + + % --------------------------------------- Test accumulate_sqw + % --------------------------------------- + + % Create some sqw files against which to compare the output of + % accumulate_sqw + % --------------------------------------------------------------------------- + [~,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs]=unpack(obj); + + [tmp_files,grid1,pix_range14]=gen_sqw (obj.spe_file([1,4]), '', sqw_file_14, efix([1,4]),... + emode, alatt, angdeg, u, v, psi([1,4]), omega([1,4]), dpsi([1,4]), gl([1,4]), gs([1,4])); + for i=1:numel(tmp_files) + assertFalse(is_file(tmp_files{i})); + end + + % Now use accumulate sqw ---------------------- + obj.proj.u=u; + obj.proj.v=v; + + spe_accum={obj.spe_file{1},'','',obj.spe_file{4}}; + present = cellfun(@(x)~isempty(x),spe_accum,'UniformOutput',true); + fin = spe_accum(present); + function ff = file_rename(ff) + [fp,fn,~]=fileparts(ff); + ff=fullfile(fp,[fn,'.tmp']); + end + tmp_fls = cellfun(@file_rename,fin,'UniformOutput',false); + clobT = onCleanup(@()obj.delete_files(tmp_fls)); + + [~,grid2,acc_pix_range14]=accumulate_sqw (spe_accum, '', sqw_file_accum,efix(1:4), ... + emode, alatt, angdeg, u, v, psi(1:4), omega(1:4), dpsi(1:4), gl(1:4), gs(1:4),'clean'); + assertEqual(grid1,grid2); + + if not(obj.save_output) + assertElementsAlmostEqual(pix_range14,acc_pix_range14) + end + % the files are different and binned into different ranges, but + % actual cuts are the same as pixels included in cuts are the + % same + [ok,mess,w2_14]=is_cut_equal(sqw_file_14,sqw_file_accum,... + obj.proj,[-1.5,0.025,0],[-2.1,-1.9],[-0.5,0.5],[-inf,inf]); + assertTrue(ok,['Cuts from gen_sqw output and accumulate_sqw are not the same',mess]); + + % Test against saved or store to save later + obj.assertEqualToTolWithSave(w2_14,'ignore_str',true,'tol',1.e-7); + + + end + % + function test_accumulate_and_combine1to4(obj,varargin) + if obj.skip_test + skipTest(fprintf('test_accumulate_and_combine1to4_%s is disabled',obj.test_pref)); + end + if nargin> 1 + % running in single test method mode. + obj.setUp(); + co1 = onCleanup(@()obj.tearDown()); + end + %------------------------------------------------------------- + % this test works with 1 tmp file so should not be run in + % parpool mode. + hpc = hpc_config; + comb_code = hpc.combine_sqw_using; + if strcmpi(comb_code,'mpi_code') + return; + end + + % build test files if they have not been build + obj=build_test_files(obj); + file_pref = obj.test_pref; + sqw_file_accum=fullfile(obj.working_dir,['test_accum_and_comb14',file_pref,'.sqw']); % output sqw file + if ~obj.save_output + co2=onCleanup(@()obj.delete_files(sqw_file_accum)); + end + + spe_names = obj.spe_file([1,4,5,6]); + for i=1:numel(spe_names) + [fp,fn,~] = fileparts(spe_names{i}); + if exist(fullfile(fp,[fn,'.tmp']),'file') == 2 + delete(fullfile(fp,[fn,'.tmp'])); + end + end + + new_names = gen_sqw_accumulate_sqw_tests_common.rename_file_list(spe_names(3:4),'.tnxs'); + co3 = onCleanup(@()gen_sqw_accumulate_sqw_tests_common.rename_file_list(new_names,'.nxspe')); + + % Allow nan-s and inf-s to keep masked pixels, to ensure + % consistency of estimated and actual pix_ranges + co5 = set_temporary_config_options(hor_config, ... + 'ignore_nan', false, ... + 'ignore_inf', false ... + ); + + % --------------------------------------- Test accumulate_sqw + % --------------------------------------- + + % Create some sqw files against which to compare the output of + % accumulate_sqw + % --------------------------------------------------------------------------- + [~,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs]=unpack(obj,[1,4,5,6]); + + % Now use accumulate sqw ---------------------- + % This is the case when the pixel_range is incorrect until the + % whole run is completed and final pixel_range equal to the + % img_db_range + [~,~,pix_range_f14]=accumulate_sqw(spe_names, '', sqw_file_accum, ... + efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs); + % new file have been reduced. Add it to the file list + gen_sqw_accumulate_sqw_tests_common.rename_file_list(new_names{1},'.nxspe'); + ldr = sqw_formats_factory.instance().get_loader(sqw_file_accum); + dat = ldr.get_data(); + clear ldr; + img_db_range1 = dat.img_range; + + % add new file to the list of the files + [~,~,pix_range_f145]=accumulate_sqw(spe_names, '', sqw_file_accum, ... + efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs); + assertTrue(all(pix_range_f14(1,:)>=pix_range_f145(1,:))); + assertTrue(all(pix_range_f14(2,:)<=pix_range_f145(2,:))); + + ldr = sqw_formats_factory.instance().get_loader(sqw_file_accum); + dat = ldr.get_data(); + clear ldr; + img_db_range2 = dat.img_range; + assertEqual(img_db_range1,img_db_range2); + + + gen_sqw_accumulate_sqw_tests_common.rename_file_list(new_names{2},'.nxspe'); + [tmp_fls,~,pix_range_f1456]=accumulate_sqw(spe_names, '', sqw_file_accum, ... + efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs); + assertTrue(all(pix_range_f145(1,:)>=pix_range_f1456(1,:))); + assertTrue(all(pix_range_f145(2,:)<=pix_range_f1456(2,:))); + % + clAccFls=onCleanup(@()obj.delete_files(tmp_fls)); + % + ldr = sqw_formats_factory.instance().get_loader(sqw_file_accum); + dat = ldr.get_data(); + clear ldr; + img_db_range3 = dat.img_range; + assertEqual(img_db_range1,img_db_range3); + + % img_db_range wider then pix_range because of energies estimate + % for initially missed runfiles from existing runfiles is wider + % then actual energy range for the existing runfiles. + assertTrue(all(pix_range_f1456(1,1:4)>=img_db_range3(1,:))); + assertTrue(all(pix_range_f1456(2,1:4)<=img_db_range3(2,:))); + + %---------------------------- + c_proj = struct('u',u,'v',v); + w2_1456=cut_sqw(sqw_file_accum,c_proj ,[-1.5,0.025,0],[-2.1,-1.9],[-0.5,0.5],[-Inf,Inf]); + + % Test against saved or store to save later + obj.assertEqualToTolWithSave(w2_1456,'ignore_str',true,'tol',1.e-7); + end + + function test_accumulate_sqw1456(obj,varargin) + %------------------------------------------------------------- + if obj.skip_test + skipTest(fprintf('test_accumulate_sqw1456_%s is disabled',obj.test_pref)); + end + if nargin> 1 + % running in single test method mode. + obj.setUp(); + co1 = onCleanup(@()obj.tearDown()); + end + %------------------------------------------------------------- + % Do not delete tmp file for future accumulation + + co2 = set_temporary_config_options(hor_config, 'delete_tmp', false); + + %------------------------------------------------------------- + + % build test files if they have not been build + obj=build_test_files(obj); + file_pref = obj.test_pref; + sqw_file_accum=fullfile(obj.working_dir,['test_sqw_accum_sqw1456_',file_pref, '.sqw']); + sqw_file_1456=fullfile(obj.working_dir,['test_sqw_1456_',file_pref, '.sqw']); + + if ~obj.save_output + cleanup_obj1=onCleanup(@()obj.delete_files(sqw_file_1456,sqw_file_accum)); + end + % --------------------------------------- Test accumulate_sqw + % --------------------------------------- + + % Create some sqw files against which to compare the output of + % accumulate_sqw + % --------------------------------------------------------------------------- + [~,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs]=unpack(obj); + hc = hpc_config; + disp(hc); + + [tmp_fls,~,pix_range1456]=... + gen_sqw (obj.spe_file([1,4,5,6]), '',sqw_file_1456, efix([1,4,5,6]),... + emode, alatt, angdeg, u, v,... + psi([1,4,5,6]), omega([1,4,5,6]), dpsi([1,4,5,6]), gl([1,4,5,6]), gs([1,4,5,6])); + for i=1:numel(tmp_fls) + assertTrue(is_file(tmp_fls{i})); + end + ldr = sqw_formats_factory.instance().get_loader(sqw_file_1456); + dat = ldr.get_data(); + data_range = ldr.get_data_range(); + clear ldr; + img_db_range1 = dat.img_range; + assertElementsAlmostEqual(data_range,pix_range1456,'relative',4*eps('single')); + assertElementsAlmostEqual(data_range(:,1:4),img_db_range1,'relative',1.e-4); + + + % Now use accumulate sqw ---------------------- + + spe_accum={obj.spe_file{1},'','',obj.spe_file{4},obj.spe_file{5},obj.spe_file{6}}; + [tmp_fls_t,~,acc_pix_range1456,all_sqw]=... + accumulate_sqw (spe_accum, '', sqw_file_accum,efix,... + emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs); + clFls = onCleanup(@()obj.delete_files(tmp_fls_t)); + for i=1:numel(tmp_fls) + assertTrue(is_file(tmp_fls_t{i})); + end + + assertElementsAlmostEqual(all_sqw.pix.data_range,acc_pix_range1456,'relative',4*eps('single')); + + img_db_range2 = all_sqw.data.img_range; + % img_db_range in second case is wider then in the first, as + % additional ranges were added from missing files + assertTrue(all(img_db_range1(1,:)>=img_db_range2(1,:))); + assertTrue(all(img_db_range1(2,:)<=img_db_range2(2,:))); + % but pixel ranges are actual pixel ranges which are the same + assertEqualToTol(pix_range1456,acc_pix_range1456,[2.e-7,2.e-7]); + + proj_o = struct('u',u,'v',v); + [ok,mess,w2_1456]=is_cut_equal(sqw_file_1456,sqw_file_accum,... + proj_o,[-1.5,0.025,0],[-2.1,-1.9],[-0.5,0.5],[-Inf,Inf]); + assertTrue(ok,['Cuts from gen_sqw output and accumulate_sqw are not the same: ',mess]) + % + clear all_sqw; + clear co2; + % Test against saved or store to save later + obj.assertEqualToTolWithSave(w2_1456,'ignore_str',true,'tol',1.e-7); + + end + % + function test_accumulate_sqw11456(obj,varargin) + %------------------------------------------------------------- + if obj.skip_test + skipTest(fprintf('test_accumulate_sqw11456_%s is disabled',obj.test_pref)); + end + if nargin> 1 + % running in single test method mode. + obj.setUp(); + co1 = onCleanup(@()obj.tearDown()); + end + + %------------------------------------------------------------- + % Do not delete tmp file for future accumulation + + co2 = set_temporary_config_options(hor_config, 'delete_tmp', false); + + %------------------------------------------------------------- + + % build test files if they have not been build + obj=build_test_files(obj); + file_pref = obj.test_pref; + sqw_file_accum=fullfile(obj.working_dir,['test_sqw_acc_sqw11456_',file_pref, '.sqw']); + sqw_file_11456=fullfile(obj.working_dir,['test_sqw_11456_',file_pref, '.sqw']); + cleanup_obj1=onCleanup(@()obj.delete_files(sqw_file_11456,sqw_file_accum)); + + % --------------------------------------- Test accumulate_sqw + % --------------------------------------- + % Create some sqw files against which to compare the output of + % accumulate_sqw + % --------------------------------------------------------------------------- + [~,efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs]=unpack(obj); + spe_selected = obj.spe_file([1,1,4,5,6]); + + + [tmp_files,grid_size1,pix_range1]=gen_sqw (spe_selected,... + '', sqw_file_11456, efix([1,3,4,5,6]), ... + emode, alatt, angdeg, u, v, psi([1,3,4,5,6]),... + omega([1,3,4,5,6]), dpsi([1,3,4,5,6]), gl([1,3,4,5,6]), gs([1,3,4,5,6]),... + 'replicate'); + clTmpDel = onCleanup(@()obj.delete_files(tmp_files)); + assertEqual(exist(sqw_file_11456,'file'),2) + for i=1:numel(tmp_files) + % files left because delete_tmp is true + assertTrue(is_file(tmp_files{i})); + end + + clear co2; + co2 = set_temporary_config_options(hor_config, 'delete_tmp', true,'log_level',1); + clWarn = set_temporary_warning('off','HORACE:push_warning','HORACE:valid_tmp_files_exist'); + % Now use accumulate sqw ---------------------- + obj.proj.u=u; + obj.proj.v=v; + + % Repeat a file with 'replicate' + spe_accum={obj.spe_file{1},'',obj.spe_file{1},obj.spe_file{4},... + obj.spe_file{5},obj.spe_file{6}}; + warning('HORACE:push_warning','push warning issued to ensure correct warning will appear below'); + [tmp_fls,~,pix_range2]=accumulate_sqw (spe_accum, '',... + sqw_file_accum,efix, emode, alatt, angdeg, u, v, psi,... + omega, dpsi, gl, gs,... + 'replicate'); %grid_size1,pix_range1, + % these files have different pixel coordinate ranges so + % existing tmp will not be reused + [~,warn_id] = lastwarn; + % Check on no warning on reusing existing files as this warning + % should not be issued in accumulate mode + files are not reused + % anyway + assertFalse(isequal(warn_id,'HORACE:valid_tmp_files_exist')) + + % Not all provided spe files were present so tmp are not deleted + % after sqw file was generated + for i=1:numel(tmp_files) + assertTrue(is_file(tmp_fls{i})); + end + + + assertEqualToTol(pix_range1,pix_range2,[2.e-7,2.e-7]); + % Estimated ranges the pixels rebinned onto are different but + % actual pix ranges are the same as the data files are the same. + % This is why the cut have to be made within the specified + % ranges to get the same result + [ok,mess,w2_11456,w2_11456acc]=... + is_cut_equal(sqw_file_11456,sqw_file_accum,obj.proj,... + [-1.5,0.025,0],[-2.1,-1.9],[-0.5,0.5],[2,35]); + if ~ok + acolor('b'); + plot(w2_11456); + acolor('r'); + pd(w2_11456acc); + keep_figure; + end + assertTrue(ok,['Cuts from gen_sqw output and accumulate_sqw are not the same',mess]); + % Test against saved or store to save later + obj.assertEqualToTolWithSave(w2_11456,'ignore_str',true,'tol',1.e-7); + + + if obj.save_output + return; + end + % Accumulate nothing, all files already accumulated. No warning + % as this is accumulate_sqw mode + spe_accum={obj.spe_file{1},'',obj.spe_file{1},obj.spe_file{4},obj.spe_file{5},obj.spe_file{6}}; + [~,grid_size,pix_range]=accumulate_sqw (spe_accum, '', sqw_file_accum,... + efix, emode, alatt, angdeg, u, v, psi, omega, dpsi, gl, gs, 'replicate'); + + assertElementsAlmostEqual(pix_range,pix_range2,'relative',4*eps('single')); + assertEqual(grid_size,grid_size1); + + [ok,mess]=is_cut_equal(sqw_file_11456,sqw_file_accum,obj.proj,... + [-1.5,0.025,0],[-2.1,-1.9],[-0.5,0.5],[2,35]); + + assertTrue(ok,['Cuts from gen_sqw output and accumulate_sqw are not the same: ',... + mess]); + clear cleanup_obj1; + clear clTmpDel; + end + % + end +end diff --git a/herbert_core/classes/data_loaders/@loader_nxspe/loader_nxspe.m b/herbert_core/classes/data_loaders/@loader_nxspe/loader_nxspe.m index e652696390..73cf87fd70 100644 --- a/herbert_core/classes/data_loaders/@loader_nxspe/loader_nxspe.m +++ b/herbert_core/classes/data_loaders/@loader_nxspe/loader_nxspe.m @@ -141,7 +141,9 @@ obj.nexus_instrument_ = obj.read_instrument_info_(); catch ME if strcmp(ME.identifier, 'HERBERT:loader_nxspe:missing_instrument_fields') - warning(ME.identifier, ME.message); + warning(ME.identifier, '%s', ME.message); + else + rethrow(ME); end % Ignore all other errors; instrument info not guaranteed to % be in all nxspe files; its absence is not an error. @@ -303,6 +305,7 @@ end function moderator = read_inst_moderator_(obj, ds) % Construct an IX_moderator from a NeXus data structure + moderator_described = true; if isfield(ds.moderator, 'pulse_shape') pulse_model = 'table'; parameters = {ds.moderator.pulse_shape.Time.value ... @@ -311,12 +314,17 @@ pulse_model = ds.moderator.empirical_pulse_shape.type.value; parameters = ds.moderator.empirical_pulse_shape.data.value; else - error('HERBERT:loader_nxspe:invalid_moderator', ... + moderator_described = false; + warning('HERBERT:loader_nxspe:invalid_moderator', ... 'moderator model in instrument info not understandable by Horace.'); end - moderator = IX_moderator(abs(ds.moderator.transforms.MOD_T_AXIS.value), ... - ds.moderator.transforms.MOD_R_AXIS.value, ... - pulse_model, parameters); + if moderator_described + moderator = IX_moderator(abs(ds.moderator.transforms.MOD_T_AXIS.value), ... + ds.moderator.transforms.MOD_R_AXIS.value, ... + pulse_model, parameters); + else + moderator = IX_inst.fill_missing_moderator(ds); + end end function instrument = read_fermi_inst_(obj, ds, src, mod) % Construct an IX_inst_DGfermi from a NeXus data structure diff --git a/herbert_core/classes/instrument_classes/instrument/@IX_fermi_chopper/private/get_pulse_props_.m b/herbert_core/classes/instrument_classes/instrument/@IX_fermi_chopper/private/get_pulse_props_.m index d83641b193..026c15d78a 100644 --- a/herbert_core/classes/instrument_classes/instrument/@IX_fermi_chopper/private/get_pulse_props_.m +++ b/herbert_core/classes/instrument_classes/instrument/@IX_fermi_chopper/private/get_pulse_props_.m @@ -37,10 +37,14 @@ else % One or more chopper parameters have been set vi = 1e6 * sqrt(ei) / obj.c_e_to_t_; % incident velocity (m/s) - - omega=2*pi*obj.frequency_; + + % attempt to fix problem with obj.frequency coming in from nxspe files populated + % with instrument information as int32. + % Cast to double for omega effective and makes gam for pk_fwhh==0 to be Inf not max(int32). + % Ordering change for pk_fwhh left in although it doesn't persuade the calc to widen to double. + omega=double(2*pi*obj.frequency_); s=2*omega*obj.curvature_; - pk_fwhh=obj.slit_width_/(2*obj.radius_*omega); + pk_fwhh=obj.slit_width_/(2*omega*obj.radius_); if phase gam=(2*obj.radius_/pk_fwhh)*abs(1/s-1./vi); diff --git a/herbert_core/classes/instrument_classes/instrument/instruments/@IX_inst/IX_inst.m b/herbert_core/classes/instrument_classes/instrument/instruments/@IX_inst/IX_inst.m index 3caaed7929..9bdfd480eb 100644 --- a/herbert_core/classes/instrument_classes/instrument/instruments/@IX_inst/IX_inst.m +++ b/herbert_core/classes/instrument_classes/instrument/instruments/@IX_inst/IX_inst.m @@ -37,6 +37,27 @@ function rv = xxx(obj) rv = isa(obj,'IX_inst'); end + + function mod = fill_missing_moderator(ds) + % temporary fix to moderator input while there is an upstream defect in moderator + % population. This code fills in "a" valid moderator. This is not intended to be used, + % rather it enables gen_sqw to reach the point where the user can populate the moderator + % manually to fix the problem. + % Currently only a default MAPS moderator is coded, other instrument defaults will be + % added as and when they are found to be needed. + if isfield(ds,'name') + if (strcmp(ds.name.value,'MAPS')) + instr = maps_instrument(300,250,'S'); + mod = instr.moderator; + else + error('HORACE:IX_inst:invalid_argument', ... + 'instruments other than MAPS not yet included'); + end + else + error('HORACE:IX_inst:invalid_argument', ... + 'name not given for instrument'); + end + end end methods %------------------------------------------------------------------ diff --git a/herbert_core/classes/instrument_classes/instrument/instruments/@IX_inst_DGfermi/IX_inst_DGfermi.m b/herbert_core/classes/instrument_classes/instrument/instruments/@IX_inst_DGfermi/IX_inst_DGfermi.m index ab964640fc..f83715096a 100644 --- a/herbert_core/classes/instrument_classes/instrument/instruments/@IX_inst_DGfermi/IX_inst_DGfermi.m +++ b/herbert_core/classes/instrument_classes/instrument/instruments/@IX_inst_DGfermi/IX_inst_DGfermi.m @@ -68,6 +68,8 @@ end function obj=set.energy(obj,val) + % Temporary fix while some variables are coming in from Matlab populated as int32 + val = double(val); obj.moderator_.energy = val; obj.fermi_chopper_.energy = val; end diff --git a/herbert_core/utilities/hdf_nexus/read_nexus_groups_recursive.m b/herbert_core/utilities/hdf_nexus/read_nexus_groups_recursive.m index adffa20db8..3a6a62978b 100644 --- a/herbert_core/utilities/hdf_nexus/read_nexus_groups_recursive.m +++ b/herbert_core/utilities/hdf_nexus/read_nexus_groups_recursive.m @@ -13,15 +13,23 @@ % switch class(filename) case 'struct' + try if ~isfield(filename, 'Filename') error('HERBERT:hdf_nexus:read_nexus_groups_recursive', 'filename and nexus_path must be strings'); elseif exist('nexus_path', 'var') error('HERBERT:hdf_nexus:read_nexus_groups_recursive', 'filename provided as struct, but nexus_path also provided'); end - + catch ME + disp(' ') + end + try dinfo = filename; filename = dinfo.Filename; + catch ME + disp(' ') + end case {'string', 'char'} + try if ~exist('nexus_path', 'var') error('HERBERT:hdf_nexus:read_nexus_groups_recursive', 'filename provided as string, but nexus_path not provided'); end @@ -31,6 +39,9 @@ else dinfo = nexus_path; end + catch ME + disp(' ') + end otherwise error('HERBERT:hdf_nexus:read_nexus_groups_recursive', ... @@ -38,14 +49,53 @@ class(filename), class(nexus_path)); end + try datastruct = read_nexus_datasets(filename, dinfo); + catch ME + disp(' ') + end for ii = 1:numel(dinfo.Groups) + try pathfields = split(dinfo.Groups(ii).Name, '/'); - datastruct.(pathfields{end}) = read_nexus_groups_recursive(filename, dinfo.Groups(ii)); + catch ME + disp(' ') + end + % while spurious fields starting 'rep_' are present in nxspe files, ignore them. + % the formation of the datastruct field has been split up to ensure that the 'rep_' + % field does not generate an invalid field error (unclear why this was happening + % but the refactor here fixes it.) + try + pfe = pathfields(end); + pfe{1} + if strncmp(pfe{1},'rep_',4) + disp(' ') + isrep = true; + elseif strcmp(pfe{1},'moderator') + disp(' ') + isrep = false; + else + isrep = false; + end + catch ME + disp(' '); + end + if ~isrep + try + substruct = read_nexus_groups_recursive(filename, dinfo.Groups(ii)); + catch ME + disp(' '); + end + try + datastruct.(pfe{1}) = substruct; + catch ME + disp(' '); + end + end end end function datasets = read_nexus_datasets(filename, nexus_struct) +try datasets = struct(); for ii = 1:numel(nexus_struct.Datasets) nexus_path = [nexus_struct.Name '/' nexus_struct.Datasets(ii).Name]; @@ -57,4 +107,7 @@ fieldname = replace(nexus_struct.Datasets(ii).Name, ' ', '_'); datasets.(fieldname) = fieldstruct; end +catch ME + disp(' ') +end end