diff --git a/API/parseClassToStructs.m b/API/parseClassToStructs.m index 4d8386272..faeb5b766 100644 --- a/API/parseClassToStructs.m +++ b/API/parseClassToStructs.m @@ -90,34 +90,16 @@ % Each cell is {1 x Inf char} -% First parse the class to a structure variable. +%% First parse the class to a structure variable. inputStruct = project.toStruct(); -%% Start by removing the cell arrays -contrastLayers = inputStruct.contrastLayers; -layerDetails = inputStruct.layerDetails; - -% If any of the contrastLayers are empty, replace the empty cells by zero -% thickness layers -for i = 1:length(contrastLayers) - thisLayer = contrastLayers{i}; - if isempty(thisLayer) - contrastLayers{i} = 0; - end -end - -% Do the same for layerDetails -if isempty(layerDetails) - layerDetails = {0}; -end - -% Pull out all the cell arrays (except priors) into one array +%% Pull out all the cell arrays (except priors) into one array problemCells{1} = inputStruct.contrastRepeatSLDs; problemCells{2} = inputStruct.data; problemCells{3} = inputStruct.dataLimits; problemCells{4} = inputStruct.simLimits; -problemCells{5} = contrastLayers; -problemCells{6} = layerDetails; +problemCells{5} = inputStruct.contrastLayers; +problemCells{6} = inputStruct.layerDetails; problemCells{7} = inputStruct.paramNames; problemCells{8} = inputStruct.backgroundParamNames; problemCells{9} = inputStruct.scalefactorNames; @@ -137,39 +119,15 @@ % Now deal with domains cell arrays if isa(project, 'domainsClass') && isa(project.domainContrasts, 'domainContrastsClass') - - domainContrastLayers = inputStruct.domainContrastLayers; - - % If any of the domainContrastLayers are empty, replace the empty - % cells by zero thickness layers - for i = 1:length(domainContrastLayers) - thisLayer = domainContrastLayers{i}; - if isempty(thisLayer) - domainContrastLayers{i} = 0; - end - end problemCells{18} = inputStruct.domainContrastRepeatSLDs; - problemCells{19} = domainContrastLayers; + problemCells{19} = inputStruct.domainContrastLayers; end if isa(project, 'domainsClass') problemCells{20} = inputStruct.domainRatioNames; end -% Fix for cell array bug with custom layers - is this needed still?? -if strcmpi(inputStruct.modelType,'custom layers') || strcmpi(inputStruct.modelType,'custom xy') - for i = 1:length(problemCells{5}) - problemCells{5}{i} = 0; - end - for i = 1:length(problemCells{19}) - problemCells{19}{i} = 0; - end - - problemCells{6} = {0}; - -end - % Also the custom files array.. if isempty(problemCells{14}) problemCells{14} = {''}; diff --git a/API/projectClass/contrastsClass.m b/API/projectClass/contrastsClass.m index 844977092..f47a4dd59 100644 --- a/API/projectClass/contrastsClass.m +++ b/API/projectClass/contrastsClass.m @@ -108,7 +108,7 @@ contrastLayers{i} = thisArray; contrastCustomFile(i) = NaN; otherwise - contrastLayers{i} = {}; + contrastLayers{i} = []; whichFile = thisContrast.model; thisContrastFileNum = find(strcmpi(whichFile, allowedNames.customFileNames)); contrastCustomFile(i) = thisContrastFileNum; diff --git a/compile/fullCompile/makeCompileArgsFull.m b/compile/fullCompile/makeCompileArgsFull.m index 07a8e13b4..f8e58ac5c 100644 --- a/compile/fullCompile/makeCompileArgsFull.m +++ b/compile/fullCompile/makeCompileArgsFull.m @@ -54,7 +54,7 @@ ARG = coder.typeof(0,[1 maxArraySize],[1 1]); ARGS_1_2{5} = coder.typeof({ARG}, [1 maxArraySize],[0 1]); ARG = coder.typeof(0,[1 10],[1 1]); -ARGS_1_2{6} = coder.typeof({ARG}, [maxArraySize 1],[1 0]); +ARGS_1_2{6} = coder.typeof({ARG}, [maxArraySize 1],[1 1]); ARG = coder.typeof('X',[1 maxArraySize],[0 1]); ARGS_1_2{7} = coder.typeof({ARG}, [1 maxArraySize],[0 1]); ARG = coder.typeof('X',[1 maxArraySize],[0 1]); diff --git a/compile/reflectivityCalculation/makeCompileArgs.m b/compile/reflectivityCalculation/makeCompileArgs.m index 1709d3d7f..a049b0159 100644 --- a/compile/reflectivityCalculation/makeCompileArgs.m +++ b/compile/reflectivityCalculation/makeCompileArgs.m @@ -54,7 +54,7 @@ ARG = coder.typeof(0,[1 maxArraySize],[1 1]); ARGS_1_2{5} = coder.typeof({ARG}, [1 maxArraySize],[0 1]); ARG = coder.typeof(0,[1 10],[1 1]); -ARGS_1_2{6} = coder.typeof({ARG}, [maxArraySize 1],[1 0]); +ARGS_1_2{6} = coder.typeof({ARG}, [maxArraySize 1],[1 1]); ARG = coder.typeof('X',[1 maxArraySize],[0 1]); ARGS_1_2{7} = coder.typeof({ARG}, [1 maxArraySize],[0 1]); ARG = coder.typeof('X',[1 maxArraySize],[0 1]); diff --git a/minimisers/NS/plotLivePointsWithEllipses.m b/minimisers/NS/plotLivePointsWithEllipses.m deleted file mode 100644 index 18b4f102e..000000000 --- a/minimisers/NS/plotLivePointsWithEllipses.m +++ /dev/null @@ -1,30 +0,0 @@ -function plotLivePointsWithEllipses(livepoints, Bs, mus) -% -% plot 2-d livepoints in unit square and bounding ellipses -% -% NOTE: mainly used for debug purposes -% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -close all - -% extract number of ellipsoids and dimension -K = size(mus, 1); -ndims = size(mus, 2); - -% plot live points and bounding ellipse for ndims=2 only -if ndims==2 - for k = 1:K - [u, v, a, b, vol] = ellipsoid_2d(squeeze(Bs(k,:,:)), squeeze(mus(k,:))); - plot(u,v,'r') - hold on - end - - plot(livepoints(:,1), livepoints(:,2), '+') - xlim([-0.5 1.5]) - ylim([-0.5 1.5]) - drawnow -end - -return - diff --git a/targetFunctions/common/groupLayers/allocateLayersForContrast.m b/targetFunctions/common/groupLayers/allocateLayersForContrast.m index aa21277f8..7dc15ee7b 100644 --- a/targetFunctions/common/groupLayers/allocateLayersForContrast.m +++ b/targetFunctions/common/groupLayers/allocateLayersForContrast.m @@ -1,4 +1,4 @@ -function thisContrastLayers = allocateLayersForContrast(contrastLayers,outParameterisedLayers,useImaginary) +function thisContrastLayers = allocateLayersForContrast(contrastLayers,outParameterisedLayers,useImaginary) % Decide which layers are needed for a particular contrast. % This function takes the master array of all layers % and extracts which parameters are required for @@ -8,30 +8,14 @@ % outParameterisedLayers - List of all the available layers % thisContrastLayers - Array detailing which layers are required for this contrast -coder.varsize('thisContrastLayers',[1000 6],[1 1]); - if useImaginary thisContrastLayers = zeros(length(contrastLayers),6); else thisContrastLayers = zeros(length(contrastLayers),5); end - for i = 1:length(contrastLayers) - if (contrastLayers(i) ~= 0) - thisLayer = outParameterisedLayers{contrastLayers(i)}; - -% % Check the length of thisLayer. If it's 6, then we have an -% % imaginary SLD defined. Combile them into one complex number at -% % this point. -% if length(thisLayer) == 6 -% compSLD = complex(thisLayer(2),thisLayer(3)); -% thisLayer = [thisLayer(1) compSLD thisLayer(4:end)]; -% end - - thisContrastLayers(i,:) = thisLayer; - else - thisContrastLayers(1,:) = []; - end +for i = 1:length(contrastLayers) + thisContrastLayers(i,:) = outParameterisedLayers{contrastLayers(i)}; end -end \ No newline at end of file +end diff --git a/targetFunctions/common/groupLayers/allocateLayersForDomainsContrast.m b/targetFunctions/common/groupLayers/allocateLayersForDomainsContrast.m deleted file mode 100644 index 52f4d9750..000000000 --- a/targetFunctions/common/groupLayers/allocateLayersForDomainsContrast.m +++ /dev/null @@ -1,28 +0,0 @@ -function thisContrastLayers = allocateLayersForDomainsContrast(contrastLayers,domainContrastLayers,outParameterisedLayers,useImaginary) -% Decide which layers are needed for a particular contrast. -% This function takes the master array of all layers -% and extracts which parameters are required for -% a particular contrast. -% -% INPUTS: -% outParameterisedLayers - List of all the available layers -% thisContrastLayers - Array detailing which layers are required for this contrast - -coder.varsize('thisContrastLayers',[1000 6],[1 1]); - -if useImaginary - thisContrastLayers = zeros(length(contrastLayers),6); -else - thisContrastLayers = zeros(length(contrastLayers),5); -end - - for i = 1:length(contrastLayers) - if (contrastLayers(i) ~= 0) - thisLayer = outParameterisedLayers{contrastLayers(i)}; - thisContrastLayers(i,:) = thisLayer; - else - thisContrastLayers(1,:) = []; - end -end - -end \ No newline at end of file diff --git a/targetFunctions/common/resampleLayers/SLDFunction.m b/targetFunctions/common/resampleLayers/SLDFunction.m index b52f9065b..3cd5974bb 100644 --- a/targetFunctions/common/resampleLayers/SLDFunction.m +++ b/targetFunctions/common/resampleLayers/SLDFunction.m @@ -26,11 +26,9 @@ if ~isempty(where) sldVal = rho(where); else - belowVals = find(x > z); - aboveVals = find(x < z); - below = belowVals(end); - above = aboveVals(1); - + below = find(x > z, 1, 'last'); + above = find(x < z, 1); + belowY = rho(below); aboveY = rho(above); diff --git a/tests/domainsTFReflectivityCalculation/domainsCustomLayersInputs.mat b/tests/domainsTFReflectivityCalculation/domainsCustomLayersInputs.mat index 237d41d14..7f517b5c7 100644 Binary files a/tests/domainsTFReflectivityCalculation/domainsCustomLayersInputs.mat and b/tests/domainsTFReflectivityCalculation/domainsCustomLayersInputs.mat differ diff --git a/tests/domainsTFReflectivityCalculation/domainsCustomLayersOutputs.mat b/tests/domainsTFReflectivityCalculation/domainsCustomLayersOutputs.mat index 21bdc36b6..deceee41d 100644 Binary files a/tests/domainsTFReflectivityCalculation/domainsCustomLayersOutputs.mat and b/tests/domainsTFReflectivityCalculation/domainsCustomLayersOutputs.mat differ diff --git a/tests/domainsTFReflectivityCalculation/domainsCustomLayersTFParams.mat b/tests/domainsTFReflectivityCalculation/domainsCustomLayersTFParams.mat index 995faf462..95bb113a3 100644 Binary files a/tests/domainsTFReflectivityCalculation/domainsCustomLayersTFParams.mat and b/tests/domainsTFReflectivityCalculation/domainsCustomLayersTFParams.mat differ diff --git a/tests/domainsTFReflectivityCalculation/domainsCustomXYInputs.mat b/tests/domainsTFReflectivityCalculation/domainsCustomXYInputs.mat index bad160d90..bfba66eb1 100644 Binary files a/tests/domainsTFReflectivityCalculation/domainsCustomXYInputs.mat and b/tests/domainsTFReflectivityCalculation/domainsCustomXYInputs.mat differ diff --git a/tests/domainsTFReflectivityCalculation/domainsCustomXYOutputs.mat b/tests/domainsTFReflectivityCalculation/domainsCustomXYOutputs.mat index 499d72383..43f8f70c8 100644 Binary files a/tests/domainsTFReflectivityCalculation/domainsCustomXYOutputs.mat and b/tests/domainsTFReflectivityCalculation/domainsCustomXYOutputs.mat differ diff --git a/tests/domainsTFReflectivityCalculation/domainsCustomXYTFParams.mat b/tests/domainsTFReflectivityCalculation/domainsCustomXYTFParams.mat index 11983e4c1..112751057 100644 Binary files a/tests/domainsTFReflectivityCalculation/domainsCustomXYTFParams.mat and b/tests/domainsTFReflectivityCalculation/domainsCustomXYTFParams.mat differ diff --git a/tests/domainsTFReflectivityCalculation/domainsStandardLayersInputs.mat b/tests/domainsTFReflectivityCalculation/domainsStandardLayersInputs.mat index de4c117d8..0db76e342 100644 Binary files a/tests/domainsTFReflectivityCalculation/domainsStandardLayersInputs.mat and b/tests/domainsTFReflectivityCalculation/domainsStandardLayersInputs.mat differ diff --git a/tests/domainsTFReflectivityCalculation/domainsStandardLayersOutputs.mat b/tests/domainsTFReflectivityCalculation/domainsStandardLayersOutputs.mat index 20c9b79ad..850a83985 100644 Binary files a/tests/domainsTFReflectivityCalculation/domainsStandardLayersOutputs.mat and b/tests/domainsTFReflectivityCalculation/domainsStandardLayersOutputs.mat differ diff --git a/tests/domainsTFReflectivityCalculation/domainsStandardLayersTFParams.mat b/tests/domainsTFReflectivityCalculation/domainsStandardLayersTFParams.mat index c04167d33..659e366e9 100644 Binary files a/tests/domainsTFReflectivityCalculation/domainsStandardLayersTFParams.mat and b/tests/domainsTFReflectivityCalculation/domainsStandardLayersTFParams.mat differ diff --git a/tests/nonPolarisedTFReflectivityCalculation/customLayersInputs.mat b/tests/nonPolarisedTFReflectivityCalculation/customLayersInputs.mat index 60a1fd61d..4882b624d 100644 Binary files a/tests/nonPolarisedTFReflectivityCalculation/customLayersInputs.mat and b/tests/nonPolarisedTFReflectivityCalculation/customLayersInputs.mat differ diff --git a/tests/nonPolarisedTFReflectivityCalculation/customLayersOutputs.mat b/tests/nonPolarisedTFReflectivityCalculation/customLayersOutputs.mat index ce5c0826a..c60a6c8ae 100644 Binary files a/tests/nonPolarisedTFReflectivityCalculation/customLayersOutputs.mat and b/tests/nonPolarisedTFReflectivityCalculation/customLayersOutputs.mat differ diff --git a/tests/nonPolarisedTFReflectivityCalculation/customLayersTFParams.mat b/tests/nonPolarisedTFReflectivityCalculation/customLayersTFParams.mat index 682b6da3e..20b9a3913 100644 Binary files a/tests/nonPolarisedTFReflectivityCalculation/customLayersTFParams.mat and b/tests/nonPolarisedTFReflectivityCalculation/customLayersTFParams.mat differ diff --git a/tests/nonPolarisedTFReflectivityCalculation/customXYInputs.mat b/tests/nonPolarisedTFReflectivityCalculation/customXYInputs.mat index 8d17d80c2..d119599fe 100644 Binary files a/tests/nonPolarisedTFReflectivityCalculation/customXYInputs.mat and b/tests/nonPolarisedTFReflectivityCalculation/customXYInputs.mat differ diff --git a/tests/nonPolarisedTFReflectivityCalculation/customXYOutputs.mat b/tests/nonPolarisedTFReflectivityCalculation/customXYOutputs.mat index ac926f116..1a6c3abc8 100644 Binary files a/tests/nonPolarisedTFReflectivityCalculation/customXYOutputs.mat and b/tests/nonPolarisedTFReflectivityCalculation/customXYOutputs.mat differ diff --git a/tests/nonPolarisedTFReflectivityCalculation/customXYTFParams.mat b/tests/nonPolarisedTFReflectivityCalculation/customXYTFParams.mat index 9f28ce4d7..6f6a3b942 100644 Binary files a/tests/nonPolarisedTFReflectivityCalculation/customXYTFParams.mat and b/tests/nonPolarisedTFReflectivityCalculation/customXYTFParams.mat differ diff --git a/tests/nonPolarisedTFReflectivityCalculation/standardLayersInputs.mat b/tests/nonPolarisedTFReflectivityCalculation/standardLayersInputs.mat index 7ae3f374b..c664a2175 100644 Binary files a/tests/nonPolarisedTFReflectivityCalculation/standardLayersInputs.mat and b/tests/nonPolarisedTFReflectivityCalculation/standardLayersInputs.mat differ diff --git a/tests/nonPolarisedTFReflectivityCalculation/standardLayersOutputs.mat b/tests/nonPolarisedTFReflectivityCalculation/standardLayersOutputs.mat index 8b0572b53..e40b4df39 100644 Binary files a/tests/nonPolarisedTFReflectivityCalculation/standardLayersOutputs.mat and b/tests/nonPolarisedTFReflectivityCalculation/standardLayersOutputs.mat differ diff --git a/tests/nonPolarisedTFReflectivityCalculation/standardLayersTFParams.mat b/tests/nonPolarisedTFReflectivityCalculation/standardLayersTFParams.mat index d8bd5d945..584509d30 100644 Binary files a/tests/nonPolarisedTFReflectivityCalculation/standardLayersTFParams.mat and b/tests/nonPolarisedTFReflectivityCalculation/standardLayersTFParams.mat differ diff --git a/tests/nonPolarisedTFReflectivityCalculation/testReflectivityCalculations.m b/tests/nonPolarisedTFReflectivityCalculation/testReflectivityCalculations.m index e8291b133..a5f370793 100644 --- a/tests/nonPolarisedTFReflectivityCalculation/testReflectivityCalculations.m +++ b/tests/nonPolarisedTFReflectivityCalculation/testReflectivityCalculations.m @@ -119,6 +119,19 @@ function loadTFParams(testCase, TFFile) methods (Test, ParameterCombination='exhaustive') %% Test High Level RAT Routines + function testRATEmpty(testCase) + % Test that RAT runs for a project without any layers. + emptyProject = projectClass(); + emptyProject.addContrast('data', 'Simulation',... + 'background', 'Background 1',... + 'bulkIn', 'SLD Air',... + 'bulkOut', 'SLD D2O',... + 'scalefactor', 'Scalefactor 1',... + 'resolution', 'Resolution 1'); + + [~, ~] = RAT(emptyProject,testCase.controlsInput); + end + function testRAT(testCase) [projectOutput, result] = RAT(testCase.project,testCase.controlsInput); diff --git a/tests/testContrastsClass.m b/tests/testContrastsClass.m index 6a1098e18..3fcb4955b 100644 --- a/tests/testContrastsClass.m +++ b/tests/testContrastsClass.m @@ -499,7 +499,7 @@ function testToStructCustomLayers(testCase) for i=1:testCase.numContrasts testCase.exampleClass.contrasts{i}.model = {'DSPC Model'}; end - testCase.exampleStruct.contrastLayers = {{} {} {}}; + testCase.exampleStruct.contrastLayers = {[] [] []}; testCase.exampleStruct.contrastCustomFile = [1 1 1]; testCase.verifyEqual(testCase.exampleClass.toStruct(testCase.allowedNames, 'custom layers', testCase.varTable), testCase.exampleStruct); diff --git a/tests/testProjectConversion/DSPCBilayerProjectClass.mat b/tests/testProjectConversion/DSPCBilayerProjectClass.mat index a3eda8e7b..18e57cf56 100644 Binary files a/tests/testProjectConversion/DSPCBilayerProjectClass.mat and b/tests/testProjectConversion/DSPCBilayerProjectClass.mat differ diff --git a/tests/testProjectConversion/monolayerVolumeModelProjectClass.mat b/tests/testProjectConversion/monolayerVolumeModelProjectClass.mat index caaf0015e..1a30c79be 100644 Binary files a/tests/testProjectConversion/monolayerVolumeModelProjectClass.mat and b/tests/testProjectConversion/monolayerVolumeModelProjectClass.mat differ diff --git a/utilities/misc/hdrload.m b/utilities/misc/hdrload.m deleted file mode 100644 index 40b35d1e6..000000000 --- a/utilities/misc/hdrload.m +++ /dev/null @@ -1,131 +0,0 @@ -function [header, data] = hdrload(file) - - -% HDRLOAD Load data from an ASCII file containing a text header. -% [header, data] = HDRLOAD('filename.ext') reads a data file -% called 'filename.ext', which contains a text header. There -% is no default extension; any extensions must be explicitly -% supplied. -% -% The first output, HEADER, is the header information, -% returned as a text array. -% The second output, DATA, is the data matrix. This data -% matrix has the same dimensions as the data in the file, one -% row per line of ASCII data in the file. If the data is not -% regularly spaced (i.e., each line of ASCII data does not -% contain the same number of points), the data is returned as -% a column vector. -% -% Limitations: No line of the text header can begin with -% a number. Only one header and data set will be read, -% and the header must come before the data. -% -% See also LOAD, SAVE, SPCONVERT, FSCANF, FPRINTF, STR2MAT. -% See also the IOFUN directory. - - -% check number and type of arguments -if nargin < 1 - error('Function requires one input argument'); -elseif ~isstr(file) - error('Input must be a string representing a filename'); -end - - -% Open the file. If this returns a -1, we did not open the file -% successfully. -fid = fopen(file); -if fid==-1 - error('File not found or permission denied'); - end - - -% Initialize loop variables -% We store the number of lines in the header, and the maximum -% length of any one line in the header. These are used later -% in assigning the 'header' output variable. -no_lines = 0; -max_line = 0; - - -% We also store the number of columns in the data we read. This -% way we can compute the size of the output based on the number -% of columns and the total number of data points. -ncols = 0; - - -% Finally, we initialise the data to []. -data = []; - - -% Start processing. -line = fgetl(fid); -if ~isstr(line) - disp('Warning: file contains no header and no data') - end; -[data, ncols, errmsg, nxtindex] = sscanf(line, '%f'); - - -% One slight problem, pointed out by Peter vanderWal: If the -% first character of the line is 'e', then this will scan as -% 0.00e+00. We can trap this case specifically by using the -% 'next index' output: in the case of a stripped 'e' the next -% index is one, indicating zero characters read. See the help -% entry for 'sscanf' for more information on this output -% parameter. We loop through the file one line at a time until -% we find some data. After that point we stop checking for -% header information. This part of the program takes most of the -% processing time, because fgetl is relatively slow (compared to -% fscanf, which we will use later). -while isempty(data)|(nxtindex==1) - no_lines = no_lines+1; - max_line = max([max_line, length(line)]); - % Create unique variable to hold this line of text information. - % Store the last-read line in this variable. - eval(['line', num2str(no_lines), '=line;']); - line = fgetl(fid); - if ~isstr(line) - disp('Warning: file contains no data') - break - end; - [data, ncols, errmsg, nxtindex] = sscanf(line, '%f'); - end % while - - -% Now that we have read in the first line of data, we can skip -% the processing that stores header information, and just read -% in the rest of the data. -data = [data; fscanf(fid, '%f')]; -fclose(fid); - - -% Create header output from line information. The number of lines -% and the maximum line length are stored explicitly, and each -% line is stored in a unique variable using the 'eval' statement -% within the loop. Note that, if we knew a priori that the -% headers were 10 lines or less, we could use the STR2MAT -% function and save some work. First, initialise the header to an -% array of spaces. -header = setstr(' '*ones(no_lines, max_line)); -for i = 1:no_lines - varname = ['line' num2str(i)]; - % Note that we only assign this line variable to a subset of - % this row of the header array. We thus ensure that the matrix - % sizes in the assignment are equal. We also consider blank - % header lines using the following IF statement. - if eval(['length(' varname ')~=0']) - eval(['header(i, 1:length(' varname ')) = ' varname ';']); - end - end % for - - -% Resize output data, based on the number of columns (as returned -% from the sscanf of the first line of data) and the total number -% of data elements. Since the data was read in row-wise, and -% MATLAB stores data in columnwise format, we have to reverse the -% size arguments and then transpose the data. If we read in -% irregularly spaced data, then the division we are about to do -% will not work. Therefore, we will trap the error with an EVAL -% call; if the reshape fails, we will just return the data as is. -eval('data = reshape(data, ncols, length(data)/ncols)'';', ''); - diff --git a/utilities/misc/jacobianEstimate/jacobianEstimate.m b/utilities/misc/jacobianEstimate/jacobianEstimate.m deleted file mode 100644 index 4a5af3219..000000000 --- a/utilities/misc/jacobianEstimate/jacobianEstimate.m +++ /dev/null @@ -1,210 +0,0 @@ -function [jac,err] = jacobianEstimate(fun,x0) -% gradest: estimate of the Jacobian matrix of a vector valued function of n variables -% usage: [jac,err] = jacobianEstimate(fun,x0) -% -% -% arguments: (input) -% fun - (vector valued) analytical function to differentiate. -% fun must be a function of the vector or array x0. -% -% x0 - vector location at which to differentiate fun -% If x0 is an nxm array, then fun is assumed to be -% a function of n*m variables. -% -% -% arguments: (output) -% jac - array of first partial derivatives of fun. -% Assuming that x0 is a vector of length p -% and fun returns a vector of length n, then -% jac will be an array of size (n,p) -% -% err - vector of error estimates corresponding to -% each partial derivative in jac. -% -% -% Example: (nonlinear least squares) -% xdata = (0:.1:1)'; -% ydata = 1+2*exp(0.75*xdata); -% fun = @(c) ((c(1)+c(2)*exp(c(3)*xdata)) - ydata).^2; -% -% [jac,err] = jacobianEstimate(fun,[1 1 1]) -% -% jac = -% -2 -2 0 -% -2.1012 -2.3222 -0.23222 -% -2.2045 -2.6926 -0.53852 -% -2.3096 -3.1176 -0.93528 -% -2.4158 -3.6039 -1.4416 -% -2.5225 -4.1589 -2.0795 -% -2.629 -4.7904 -2.8742 -% -2.7343 -5.5063 -3.8544 -% -2.8374 -6.3147 -5.0518 -% -2.9369 -7.2237 -6.5013 -% -3.0314 -8.2403 -8.2403 -% -% err = -% 5.0134e-15 5.0134e-15 0 -% 5.0134e-15 0 2.8211e-14 -% 5.0134e-15 8.6834e-15 1.5804e-14 -% 0 7.09e-15 3.8227e-13 -% 5.0134e-15 5.0134e-15 7.5201e-15 -% 5.0134e-15 1.0027e-14 2.9233e-14 -% 5.0134e-15 0 6.0585e-13 -% 5.0134e-15 1.0027e-14 7.2673e-13 -% 5.0134e-15 1.0027e-14 3.0495e-13 -% 5.0134e-15 1.0027e-14 3.1707e-14 -% 5.0134e-15 2.0053e-14 1.4013e-12 -% -% (At [1 2 0.75], jac should be numerically zero) -% -% -% See also: derivest, gradient, gradest -% -% -% Author: John D'Errico -% e-mail: woodchips@rochester.rr.com -% Release: 1.0 -% Release date: 3/6/2007 - -% get the length of x0 for the size of jac -nx = numel(x0); - -MaxStep = 100; -StepRatio = 2; - -% get fun at the center point -f0 = fun(x0); -f0 = f0(:); -n = length(f0); -if n==0 - % empty begets empty - jac = zeros(0,nx); - err = jac; - return -end - -relativedelta = MaxStep*StepRatio .^(0:-1:-25); -nsteps = length(relativedelta); - -% total number of derivatives we will need to take -jac = zeros(n,nx); -err = jac; -for i = 1:nx - x0_i = x0(i); - if x0_i ~= 0 - delta = x0_i*relativedelta; - else - delta = relativedelta; - end - - % evaluate at each step, centered around x0_i - % difference to give a second order estimate - fdel = zeros(n,nsteps); - for j = 1:nsteps - fdif = fun(swapelement(x0,i,x0_i + delta(j))) - ... - fun(swapelement(x0,i,x0_i - delta(j))); - - fdel(:,j) = fdif(:); - end - - % these are pure second order estimates of the - % first derivative, for each trial delta. - derest = fdel.*repmat(0.5 ./ delta,n,1); - - % The error term on these estimates has a second order - % component, but also some 4th and 6th order terms in it. - % Use Romberg exrapolation to improve the estimates to - % 6th order, as well as to provide the error estimate. - - % loop here, as rombextrap coupled with the trimming - % will get complicated otherwise. - for j = 1:n - [der_romb,errest] = rombextrap(StepRatio,derest(j,:),[2 4]); - - % trim off 3 estimates at each end of the scale - nest = length(der_romb); - trim = [1:3, nest+(-2:0)]; - [der_romb,tags] = sort(der_romb); - der_romb(trim) = []; - tags(trim) = []; - - errest = errest(tags); - - % now pick the estimate with the lowest predicted error - [err(j,i),ind] = min(errest); - jac(j,i) = der_romb(ind); - end -end - -end % mainline function end - -% ======================================= -% sub-functions -% ======================================= -function vec = swapelement(vec,ind,val) -% swaps val as element ind, into the vector vec -vec(ind) = val; - -end % sub-function end - -% ============================================ -% subfunction - romberg extrapolation -% ============================================ -function [der_romb,errest] = rombextrap(StepRatio,der_init,rombexpon) -% do romberg extrapolation for each estimate -% -% StepRatio - Ratio decrease in step -% der_init - initial derivative estimates -% rombexpon - higher order terms to cancel using the romberg step -% -% der_romb - derivative estimates returned -% errest - error estimates -% amp - noise amplification factor due to the romberg step - -srinv = 1/StepRatio; - -% do nothing if no romberg terms -nexpon = length(rombexpon); -rmat = ones(nexpon+2,nexpon+1); -% two romberg terms -rmat(2,2:3) = srinv.^rombexpon; -rmat(3,2:3) = srinv.^(2*rombexpon); -rmat(4,2:3) = srinv.^(3*rombexpon); - -% qr factorization used for the extrapolation as well -% as the uncertainty estimates -[qromb,rromb] = qr(rmat,0); - -% the noise amplification is further amplified by the Romberg step. -% amp = cond(rromb); - -% this does the extrapolation to a zero step size. -ne = length(der_init); -rhs = vec2mat(der_init,nexpon+2,ne - (nexpon+2)); -rombcoefs = rromb\(qromb'*rhs); -der_romb = rombcoefs(1,:)'; - -% uncertainty estimate of derivative prediction -s = sqrt(sum((rhs - rmat*rombcoefs).^2,1)); -rinv = rromb\eye(nexpon+1); -cov1 = sum(rinv.^2,2); % 1 spare dof -errest = s'*12.7062047361747*sqrt(cov1(1)); - -end % rombextrap - - -% ============================================ -% subfunction - vec2mat -% ============================================ -function mat = vec2mat(vec,n,m) -% forms the matrix M, such that M(i,j) = vec(i+j-1) -[i,j] = ndgrid(1:n,0:m-1); -ind = i+j; -mat = vec(ind); -if n==1 - mat = mat'; -end - -end % vec2mat - -