Double parameter optimisation using optimisation toolbox

Hi,

I am wondering if someone could start me off in the right direction.

I have two parameters. P1 ranges from 30:-30 and P2 from 0.00000000000000001:1.

These two parameters can be used to generate a graph. I can then have the X,Y data for that graph.

I have a second graph that i obtain from an experiment and i want to know what values of P1 and P2 will produce that graph(or one that is similar to a fine value).

So i need to sum (experimentalYpoint - generatedYpoint)^2 across all my x points. This will give me a value of how 'similar' the two are. i then optimise P1 and P2 so that the value of my sum gets as close to 0 as it can be.

I think i can use the Optimization Toolbox for this but would someone be able to point me to exactly which part of it that would be useful.

Forry for the long winded question, wanted to try make it clear.

[Torsten]

help lsqcurvefit

[Charles Mitchell-Thurston]

'given input data xdata, and the observed output ydata, where xdata and ydata are matrices or vectors, and F (x, xdata) is a matrix-valued or vector-valued function of the same size as ydata.'

In this is my xdata P1 and ydata P2?

[Torsten]

xdata = independent variable on which your "experimentalYpoints" depend (e.g. time, length or whatever)

ydata = experimentalYpoints

x = [P1,P2] = vector of parameters to be fitted

Answer 1

[Matt J]

In this is my xdata P1 and ydata P2?

xdata, ydata are what you call experimtalX and experimentalY. The unknown parameters that you are trying to solve for go in the vector x, that is x=[P1,P2].

[Charles Mitchell-Thurston]

Im probably being thick but what do i actually put them as? what do i meake P1 and P2 equal to?

[Torsten]

P1 and P2 are the unknowns. You only need to assign initial values to them - lsqcurvefit will try to optimize them such that sum (experimentalYpoint - generatedYpoint(P1,P2))^2 becomes minimal.

[Charles Mitchell-Thurston]

Right that makes sense.

I am still a bit unsure on how to set this all up though.

So x=[15,0.5] and this can just be my starting guess as its in the middle of all possible values.

curve_fcn = @(P1P2,x) P1P2(1)*(exp(-(x-P1P2(2)).^2/1e3) + exp(-(x+P1P2(2)).^2/1e3)); Does what my curve function actually is matter?

Using Bjorn's answer as well:

res_fcn = @(P1P2,xdata,ydata,fcn) (y-fcn(P1P2,xdata));

where xdata is xexperimental and ydata is yexperimental

P1P2_lsq = lsqnonlin(@(par) res_fcn(par,xdata,ydata,curve_fcn),x);

So i would need something like

x=[15,0.5]; %standard middle guess
curve_fcn = @(P1P2,x) P1P2(1)*(exp(-(x-P1P2(2)).^2/1e3) + exp(-(x+P1P2(2)).^2/1e3));
xdata = x experimentaldata %extracted from the experimental graph.
ydata = y experimenatldata

%my code that takes in the two parameters and output xdatasimulated and
%ydatasimulated

What i dont get is how i implament my code that gets the generated data into this all nd how i set it up. I have tried reading help lsqcurvefit but i still dont understand

[Bjorn Gustavsson]

Your xdata and ydata should go into the fitting-function call. You can generate, or read in from file, those variables any way you want. Perhaps something like this:

load xdata.mat
load ydata.mat

Your function of x that depends on your two parameters P1 and P2 is "some function" that you have told us nothing more about than that and that it produces curves like the one in the figure. Provided that function exists in an m-file or that you can express it as a one-line function and make a function-handle as I've exampled you have everything you need to call fminsearch, lsqnonlin or lsqcurevfit. With the assumption that f is a function named your_fcn in a file your_fcn.m you will do this:

P1P2_guess = [15,0.5]; %standard middle guess - use descriptive variable-names, avoid plain x for this case
P2P2 = lsqcurvefit(@(par,x) your_fcn(par(1),par(2),x),P1P2_guess,xdata,ydata);

HTH

[Charles Mitchell-Thurston]

So how the parameters work is that P1 and P2 are used to generate an FID using a program called SIMPSON using a console command in matlab. I then use some plotting code to fourier transform this into a 2 graph/spectrum and then extract the x and y coordinates using some code i found on mathworks a while ago.

[Torsten]

What lsqcurvefit and all other optimizers expect is that you are able to return "generatedYpoints" for parameters P1 and P2 that the optimizer suggests.

So you will have to automate the process you describe above somehow or just make a trial-and-error- fit without any optimizer.

[Bjorn Gustavsson]

@Charles Mitchell-Thurston, you will have to explain this in more detail. What does the SIMPSON program do, what input-parameters does it take? What is FID?

If these are what I think they are then you will have to do quite a bit more programming to automate the procedure and put that into a function that can produce the required products automatically from the inputs described in, for example, your_fcn.

[Charles Mitchell-Thurston]

Yer this is what i feared, that it might be more complex. SIMPSON is a simulation tool for chemistry. Where you input a P1 and P2. It then simulates a set number of crystal structures and generates a Free induction decay. Basically a sine curve that shrinks to 0 over time. This is then fourier transfromed to give a spectrum like the one in my origional question.

For Code i currently have, i have that you can give my code an excel file containing any number of P1 and their paired P2 values and it will do all of that above one after another and you will end up with the ammount same ammouunt of spectra which you can then extract the x and y data from.

[Bjorn Gustavsson]

@Charles Mitchell-Thurston, then things might not be too horrible anyway. What you'll have to do then is to write a function something like this (very skellington-code™):

function y_spec = your_SIMPSON_interfacer(P1,P2,x)
SIMP_CALL = ['SIMPSON ',sprintf('%f',P1),' ',sprintf('%f',P2),' ','tmp-out.dat']);
[Sout1,Sout2] = system(SIMP_CALL);
y_raw = read_SIMPSON_res('tmp-out.dat');
y_spec = SIMPSON_raw2spec(y_raw,x,P1,P2);
end

Not necessarily the most elegant way, but you should get the point. Since we don't have that much info on exactly how you call SIMPSON or its results are stored or what processing you do to get the spectrum these can only be described with this rough box-level calls.

[Charles Mitchell-Thurston]

To try help show whats goin on. Below is the part of my code where you take in data. And at the end you end up with inter_y/x Which is the X and Y data for each input that you give it. It is likely ugly code but it does the job.

The input to simpson has to be in the form of a RR.in file which contains all of the parmeters. the part of this that is changed for each parameter is in the code below where it simply inserts P1 and P2 into the .in file.

%%%%take in data%%%%
H = readmatrix('p1p2values.xlsx');
%%%%creates folder for use%%%%
mkdir RRfiles

%%%%creates .in files%%%%

text = fileread('RR.in');
index = strfind(text, '0p ');
text1 = text(1:index+2);
tempstr = text(index+3:end);
index = strfind(tempstr, ' ');
index = index(2);
text2 = tempstr(index:end);

for i = 1:2 %length(H);
    aniso =H(i,2);
    asym =H(i,3);
    newtext = [text1, num2str(aniso),  'p ',  num2str(asym), text2] ;
    filename = ['RRfiles/RR', num2str(i), '.in'];
    file = fopen(filename,'w');
    fprintf(file,newtext);
    fclose(file);
end

%%%%create fid files%%%%

list=length(H);

for k = 1:2 %list
    command = sprintf('simpson RRfiles/RR%g.in', k);
    status = system(command);
end
zero = 0;
for loop = 1:1000
    for i = 1:2 %list
        QNAN_text=fileread(['RR', num2str(i), '.fid']);
        strfind(QNAN_text, 'QNAN');
        if ans>zero
            fprintf 'Fail'
            command = sprintf('simpson RRfiles/RR%g.in', i);
            status = system(command);
        end
    end
end
fprintf 'all passed'
mkdir Fidfiles

for t = 1:2 %list;
    filename = ['RR', num2str(t), '.fid'];
    movefile ( filename, 'Fidfiles');
end

%%%% create fig files%%%%
fids = uigetdir ('C:\Users\charles\Documents\MATLAB\SIMPSON');
info = dir(fullfile(fids,'*.fid'));
list = {info.name}

%%%% Parameters%%%%
fu=' Hz';
dc = 32;
ls = 0;
wt = 'Gaussian'
wp = 10
zf = 256;
bo = 1;
bp = 32;

%%%% Process%%%%
list = natsortfiles(list);
for i = 1:length(list);
%%%%takes in data%%%%
    raw = importdata( fullfile(fids, list{i}));
    data = raw.data;
    data = data(:,1);
    FID = data;
    num = raw.textdata;
    num = num{3};
    Str = num;
    Str(Str < '0' | Str > '9') = [];
    sw = sscanf(Str, '%d');
    lf=-sw/2; % temporary scale
    hf=sw/2;
    
%%%%processes and plots the figure%%%%
    figure(i);
    f = figure;
    FID=normalize(FID);
    FID=dcOffset(FID,dc);
    FID=leftShift(FID,ls);
    FID=windowFID(FID,sw,wt,wp);
    SPE=FT(FID,zf);
    % SPE=phaseSpectrum0(SPE,-90.);
    SPE=baselineCorrect(SPE,sw,bo,bp);
    FREQ=getFrequency(SPE,lf,hf);
    % [th0,th1,pv]=Phasing(FREQ,SPE);
    %SPE = real(SPE)
    plotSpectrum(gca,FREQ,SPE);
    clear FID FREQ SPE lf hf sw test Str data raw
    savefig( num2str(i) + "hydrogen" );
    close all
end

%%%%Gets data%%%%

data = pwd;
dinfo = dir(fullfile(data,'*gen.fig'));
fignames = {dinfo.name};
numfig = length(fignames);

%%%% creates cells to be filled%%%%

y = cell(numfig, 1);
z = cell(numfig, 1);
inter_y = cell(numfig, 1);
inter_z = cell(numfig, 1);

%%%%code to create variable list for protons iso values%%%%


isofm = H(:,1);
isofm = round(isofm,1)


%%%%code that gets x,y data from sims%%%%

for K = 1 : numfig
  figfile = fignames{K};
  try
    fig = openfig(figfile,'invisible');
    ax = get(fig, 'CurrentAxes');
    if ~isempty(ax)
      hline = get(ax, 'Children');
        y{K} = get(hline,'XData');
        z{K} = get(hline,'YData');
        
        inter_y{K} = y{K}(1,5:124);
        inter_z{K} = abs(z{K}(1,5:124));
    end
    close(fig);
  end
end

[Bjorn Gustavsson]

At this point someone (saddly me) should tell you: take a step back, look at this code and figure out how to make it automated such that you automatically calculate your inter_y and inter_z (and inter_x?) values corresponding to any given [P1, P2] such that you can run an automatic fitting-search. Or you could just do a brute-force grid-search if you can afford to calculate everything for a sufficiently dense P1-P2 grid. I hope that you have enough time to take the first approach since that will be a more neatly designed program that can be used and adapted multiple ways in future applications.

[Charles Mitchell-Thurston]

So it is automated in the sense that i can get it to output the inter_y and inter_z (its called z here as its technically a z value for a 3D graph that can be created later on but that doesnt matter) for any values that i put in the excel file. I just dont know how to convert this to a usable function for lsqcurvefit.

Could, in a blunt way, just surround the whole thing with function ...... end and take out the uigetdir part so it involes 0 input outside of two values from a spreadsheet/matrix?

My orional idea was to simulate all possible spectra but what is P1 ranging from 0 to 30 in incements of 0.1 and P2 going from 0.01 to 1 in increments of 0.01. Each simulation times 30s-1min.

[Bjorn Gustavsson]

@Charles Mitchell-Thurston, yes pretty much so. As long as you can wrap this into a function that takes P1 and P2 as input-arguments (possibly also an x-array to set the dimensions of your output?) and returns the corresponding y-array. Whatever goes on inside that function is comparatively completely irrelevant to the fitting-function (as long as the function doesn't require user-selections or interactions). From the outside perspective it is all the same as long as you can write a function that produces the output given the input-parameters - except the computational time might vary rather much.

[Charles Mitchell-Thurston]

Right i think i am starting to undersand.

So i need something like

function [P1,P2] = TheoreticalY(y)

Where i can just input any two values for P1 and P2 and it will give me out a set of Simulated Y values. i can have a standard set of Sim X values as they are all the same. And then the rest is 'easy'

Also thank you for all the responses (Bjorn, torsten and matt) so far and the time you have spent helping me.

[Bjorn Gustavsson]

Pretty much. Except you've swapped the input and output-parameters around, I think. Ought to be:

function y = TheoreticalY(P1,P2)

Or for a design that's simpler to generalise:

function y = TheoreticalY(pars)
% Then you extract P1 P2 (and possible extended parameters) first:
P1 = pars(1);
P2 = pars(2);
if numel(pars) > 2 % this way you can expand the number of free parameters
    % and handle such cases without too much extra work in setting up the
    % fitting
    P3 = pars(3);
else
    P3 = exp(1i*(sqrt(5)+1)/2);
end

[Charles Mitchell-Thurston]

Yer i just realised that myself as i am trying to top and tail my code with function and end. I have tested the code below and it works for single inputs of P1 and P2 now.

The output of all of this code is inter_z which is a 1x1 cell which contails my 1x120 double which is all of my y values for my plot.

For writing that actual function, function y = TheoreticalY(P1,P2) i assume looks like a start. Do i have to somehow set that 'y' is inter_z? Do i have to clear all variables produced outside of inter_z perhaps?

P1 = 17
P2 = 0.1
%function y = TheoreticalY(P1,P2)

%%%%creates folder for use%%%%
mkdir RRfiles

%%%%creates .in files%%%%

text = fileread('RR.in');
index = strfind(text, '0p ');
text1 = text(1:index+2);
tempstr = text(index+3:end);
index = strfind(tempstr, ' ');
index = index(2);
text2 = tempstr(index:end);

newtext = [text1, num2str(P1),  'p ',  num2str(P2), text2] ;
filename = ['RRfiles/RR.in'];
file = fopen(filename,'w');
fprintf(file,newtext);
fclose(file);


%%%%create fid files%%%%

command = sprintf('simpson RRfiles/RR.in');
status = system(command);
zero = 0;
for loop = 1:1000
    QNAN_text=fileread(['RR.fid']);
    strfind(QNAN_text, 'QNAN');
    if ans>zero
        fprintf 'Fail'
        command = sprintf('simpson RRfiles/RR.in');
        status = system(command);
    end
end
fprintf 'all passed'
mkdir Fidfiles

filename = ['RR.fid'];
movefile ( filename, 'Fidfiles');


%%%% create fig files%%%%
fids = 'C:\Users\charl\Documents\MATLAB stuff\CSA fitting\Fidfiles';
info = dir(fullfile(fids,'*.fid'));
list = {info.name}

%%%% Parameters%%%%
fu=' Hz';
dc = 32;
ls = 0;
wt = 'Gaussian'
wp = 10
zf = 256;
bo = 1;
bp = 32;

%%%% Process%%%%

%%%%takes in data%%%%
    raw = importdata( fullfile(fids, list{1}));
    data = raw.data;
    data = data(:,1);
    FID = data;
    num = raw.textdata;
    num = num{3};
    Str = num;
    Str(Str < '0' | Str > '9') = [];
    sw = sscanf(Str, '%d');
    lf=-sw/2; % temporary scale
    hf=sw/2;
    
%%%%processes and plots the figure%%%%
    figure(1);
    f = figure;
    FID=normalize(FID);
    FID=dcOffset(FID,dc);
    FID=leftShift(FID,ls);
    FID=windowFID(FID,sw,wt,wp);
    SPE=FT(FID,zf);
    % SPE=phaseSpectrum0(SPE,-90.);
    SPE=baselineCorrect(SPE,sw,bo,bp);
    FREQ=getFrequency(SPE,lf,hf);
    % [th0,th1,pv]=Phasing(FREQ,SPE);
    %SPE = real(SPE)
    plotSpectrum(gca,FREQ,SPE);
    clear FID FREQ SPE lf hf sw test Str data raw
    savefig( num2str(1) + "hydrogen" );
    close all


%%%%Gets data%%%%

data = pwd;
dinfo = dir(fullfile(data,'*gen.fig'));
fignames = {dinfo.name};
numfig = length(fignames);

%%%% creates cells to be filled%%%%

y = cell(numfig, 1);
z = cell(numfig, 1);
inter_y = cell(numfig, 1);
inter_z = cell(numfig, 1);

%%%%code to create variable list for protons iso values%%%%


isofm = P1
isofm = round(isofm,1)


%%%%code that gets x,y data from sims%%%%

for K = 1 : numfig
  figfile = fignames{K};
  try
    fig = openfig(figfile,'invisible');
    ax = get(fig, 'CurrentAxes');
    if ~isempty(ax)
      hline = get(ax, 'Children');
        y{K} = get(hline,'XData');
        z{K} = get(hline,'YData');
        
        inter_y{K} = y{K}(1,5:124);
        inter_z{K} = abs(z{K}(1,5:124));
    end
    close(fig);
  end
end

[Charles Mitchell-Thurston]

I realised i was still saving it as a live script not a function. I think i have now done the hardest part of my project. Again many thanks to all, there may be another ask on matlab at a later date if i get stuck again but i think the hardest part may be over (he says not knowing the future....)

function y = TheoreticalY(P1,P2)

%%%%creates folder for use%%%%
mkdir RRfiles

%%%%creates .in files%%%%

text = fileread('RR.in');
index = strfind(text, '0p ');
text1 = text(1:index+2);
tempstr = text(index+3:end);
index = strfind(tempstr, ' ');
index = index(2);
text2 = tempstr(index:end);

newtext = [text1, num2str(P1),  'p ',  num2str(P2), text2] ;
filename = ['RRfiles/RR.in'];
file = fopen(filename,'w');
fprintf(file,newtext);
fclose(file);


%%%%create fid files%%%%

command = sprintf('simpson RRfiles/RR.in');
status = system(command);
zero = 0;
for loop = 1:1000
    QNAN_text=fileread(['RR.fid']);
    strfind(QNAN_text, 'QNAN');
    if ans>zero
        fprintf 'Fail'
        command = sprintf('simpson RRfiles/RR.in');
        status = system(command);
    end
end
fprintf 'all passed'
mkdir Fidfiles

filename = ['RR.fid'];
movefile ( filename, 'Fidfiles');


%%%% create fig files%%%%
fids = 'C:\Users\charl\Documents\MATLAB stuff\CSA fitting\Fidfiles';
info = dir(fullfile(fids,'*.fid'));
list = {info.name}

%%%% Parameters%%%%
fu=' Hz';
dc = 32;
ls = 0;
wt = 'Gaussian'
wp = 10
zf = 256;
bo = 1;
bp = 32;

%%%% Process%%%%

%%%%takes in data%%%%
    raw = importdata( fullfile(fids, list{1}));
    data = raw.data;
    data = data(:,1);
    FID = data;
    num = raw.textdata;
    num = num{3};
    Str = num;
    Str(Str < '0' | Str > '9') = [];
    sw = sscanf(Str, '%d');
    lf=-sw/2; % temporary scale
    hf=sw/2;
    
%%%%processes and plots the figure%%%%
    figure(1);
    f = figure;
    FID=normalize(FID);
    FID=dcOffset(FID,dc);
    FID=leftShift(FID,ls);
    FID=windowFID(FID,sw,wt,wp);
    SPE=FT(FID,zf);
    % SPE=phaseSpectrum0(SPE,-90.);
    SPE=baselineCorrect(SPE,sw,bo,bp);
    FREQ=getFrequency(SPE,lf,hf);
    % [th0,th1,pv]=Phasing(FREQ,SPE);
    %SPE = real(SPE)
    plotSpectrum(gca,FREQ,SPE);
    clear FID FREQ SPE lf hf sw test Str data raw
    savefig( num2str(1) + "hydrogen" );
    close all


%%%%Gets data%%%%

data = pwd;
dinfo = dir(fullfile(data,'*gen.fig'));
fignames = {dinfo.name};
numfig = length(fignames);

%%%% creates cells to be filled%%%%

y = cell(numfig, 1);
z = cell(numfig, 1);
inter_y = cell(numfig, 1);
inter_z = cell(numfig, 1);

%%%%code to create variable list for protons iso values%%%%


isofm = P1
isofm = round(isofm,1)


%%%%code that gets x,y data from sims%%%%

for K = 1 : numfig
  figfile = fignames{K};
  try
    fig = openfig(figfile,'invisible');
    ax = get(fig, 'CurrentAxes');
    if ~isempty(ax)
      hline = get(ax, 'Children');
        y{K} = get(hline,'XData');
        z{K} = get(hline,'YData');
        
        inter_y{K} = y{K}(1,5:124);
        inter_z{K} = abs(z{K}(1,5:124));
    end
    close(fig);
  end
end
end

[Bjorn Gustavsson]

Great. In the function you'll have to assign the relevant result to the y-variable.

Perhaps end the function with:

y = cell2mat(inter_y{K});

[Charles Mitchell-Thurston]

Loverly

Answer 2

[Bjorn Gustavsson]

The standard fitting method I use is:

% Your curve-function definition (I just mock one up, if you have it defined in a 
% matlab-function you don't need this step)
curve_fcn = @(P1P2,x) P1P2(1)*(exp(-(x-P1P2(2)).^2/1e3) + exp(-(x+P1P2(2)).^2/1e3));
% a sum-of-squares error-function for the fitting
err_fcn = @(P1P2,x,y,fcn) sum((y-fcn(P1P2,x)).^2);
% or a residual between the observed y and your curve-function
res_fcn = @(P1P2,x,y,fcn) (y-fcn(P1P2,x));
P1P2_guess = [1,2]; % Initial guess for the parameters P1 and P2, (adjust to be good values)
% Parameter-fitting with fminsearch
P1P2_fms = fminsearch(@(par) err_fcn(par,x,y,curve_fcn),P1P2_guess);
% Parameter-fitting with lsqnonlin:
P1P2_lsq = lsqnonlin(@(par) res_fcn(par,x,y,curve_fcn),P1P2_guess);

That should give you good parameters for P1 and P2 as the first and second element of P1P2_fms and/or P1P2_lsq.

I find these functions easier to work with since this is more transparend and I know what goes in (old dogs and new tricks are said to not go well...).

HTH

[Matt J]

???

But the only difference between your workflow and lsqcurvefit, is that lsqcurvefit takes your curv_fcn directly:

curve_fcn = @(P1P2,x) P1P2(1)*(exp(-(x-P1P2(2)).^2/1e3) + exp(-(x+P1P2(2)).^2/1e3));

P1P2_guess = [1,2]; % Initial guess for the parameters P1 and P2, (adjust to be good values)

P1P2_lsq = lsqcurvefit(curve_fcn,P1P2_guess, x,y);

Moreover, in either workflow, you will inevitably need to provide curve_fcn for when it comes time to plot the fit.

[Bjorn Gustavsson]

@Matt J, you're right (again). It is probably only due to habit, I started using fmins/fminsearch (and fminsearchbnd) and lsqnonlin and never bothered starting to use lsqcurvefit because I had gotten used to the setups for the others. You teach me a lot about myself recently, I don't like what I learn but I do appreciate it...