initlalize several variables at once

2 views (last 30 days)
Benjamin
Benjamin on 19 Mar 2019
Answered: Matt J on 20 Mar 2019
I currently have a code where I am initializing many parameters, that will be solved for later.
I do this by:
C1=params(1);
C2=params(2);
C3=params(3);
C4=params(4);
C5=params(5);
C6=params(6);
C7=params(7);
C8=params(8);
C9=params(9);
C10=params(10);
C11=params(11);
C12=params(12);
C13=params(13);
There's got to be an easier way...
Additionally, I set my intial guesses like:
params0=[0.1,0.1,0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1];
Is there a way to more easily set this, assuming they all have the same initial guess? Typing them all out that way is tedious, and my code is constantly changing. How can I automate it?

Accepted Answer

Matt J
Matt J on 19 Mar 2019
Edited: Matt J on 19 Mar 2019
If all you're trying to do is polynomial surface fits of various orders, I would just use the Curve Fitting Toolbox, if you have it
f = fit([x,y],z, 'poly55')
See also,
  22 Comments
Matt J
Matt J on 20 Mar 2019
Edited: Matt J on 20 Mar 2019
I would post that as a new question. Our back and forth has already deviated considerably from your originally posted topic.
Benjamin
Benjamin on 20 Mar 2019
Yeah, I posted a new question. Thanks. I accepted the answer here as you are right, we have deviated a lot

Sign in to comment.

More Answers (3)

Matt J
Matt J on 19 Mar 2019
Edited: Matt J on 19 Mar 2019
Is there a way to more easily set this, assuming they all have the same initial guess?
One way:
params0=linspace(0.1,0.1,13)
Another is:
params0=repelem(0.1,1,13);
There's got to be an easier way...
If you have a very large vector of variables, it does not make sense to assign them to separate variables. You should probably be manipulating them as a vector. But if you must index the variables individually, I would just rename params to something shorter,
C=params;
and then refer to the variables via indexing C(1), C(2),..C(13) throughout the code. Although, it is then unclear why you chose to name the vector params in the first place...
  4 Comments
Matt J
Matt J on 19 Mar 2019
Edited: Matt J on 19 Mar 2019
Perhaps instead of zeroth, first,etc... you generate a cell array called "orders" like so:
T=[0 0; 1 0; 0 1]; %table of exponent combinations
[orders{1:max(d)+1}]=deal(0); %initialize
for k=1:numel(C)
i=T(k,1); j=T(k,2);
orders{i+j+1}=orders{i+j+1} + C(k) * X.^i .* Y.^j ;
end
Benjamin
Benjamin on 19 Mar 2019
Edited: Benjamin on 19 Mar 2019
I'm currently doing this: I'm not exactly sure how to replace it with your code. The way I currently have it allows be to check the validity of my fit starting at zeroth order, and then I can individually add in higher order. How can I use your code, which seems much more efficient, to do what I am doing here? Where I just change the order to say, 4, or 5, and accomplishes what I am below? I'm finding that a 5th order polynomial fits my data really well. Unfortunately, the surface extends beyond my data, and I am not sure how to automatically get rid of those values, in order to make my surface look nice. I know you said the accuracy calcuation is not affected, but my graph would look so much better if the surface could be chopped off somehow in regions where there is no data.
zeroth = C(1) .* (X.^0) .* (Y.^0);
first = C(2) .* (X.^1) .* (Y.^0) + C(3) * (X.^0) .* (Y.^1);
second = C(4) .* (X.^2) .* (Y.^0) + C(5) * (X.^1) .* (Y.^1) + C(6) * (X.^0) .* (Y.^2);
third = C(7) .* (X.^3) .* (Y.^0) + C(8) * (X.^2) .* (Y.^1) + C(9) * (X.^1) .* (Y.^2) + C(10) * (X.^0) .* (Y.^3);
fourth = C(11) .* (X.^4) .* (Y.^0) + C(12) * (X.^3) .* (Y.^1) + C(13) * (X.^2) .* (Y.^2) + C(14) * (X.^1) .* (Y.^3)+ C(15) * (X.^0) .* (Y.^4);
fifth = C(16) .* (X.^5) .* (Y.^0) + C(17) * (X.^4) .* (Y.^1) + C(18) * (X.^3) .* (Y.^2) + C(19) * (X.^2) .* (Y.^3)+ C(20) * (X.^1) .* (Y.^4)+ C(21) * (X.^0) .* (Y.^5);
poly_guess = zeroth+first+second+third+fourth;

Sign in to comment.


Matt J
Matt J on 20 Mar 2019
Edited: Matt J on 20 Mar 2019
One other thing I would point out is the model function you have shown us here is linear in the coefficients. This means that it can be represented as a matrix multiplication and the entire fit can be done with simple linear algebra instead of an iterative solver like lsqcurvefit.
To faciliate this you could download my func2mat utility and then do as follows
A=func2mat(@(p) modelfun(p,xdata), params0);
coefficients =A\g(:);
That's it!
  4 Comments
Matt J
Matt J on 20 Mar 2019
Nothing. Whatever modelfun() you were planning to use with lsqcurvefit, you simply use instead in the above 2 lines of code.
Benjamin
Benjamin on 20 Mar 2019
Hmmm, I can't seem to figure it out... Here is my code (in all of its glory). What do I need to do to change from lsq to matrix multiplication? Btw, the lsq method seems to be working really well. The fits are looking amazing, and I appreciate all your help with getting that to work!
clear
clc
close all
num=xlsread('C:\example.xlsx',1,'A2:F18110');
eta = num(:,3);
r = num(:,4);
g = num(:,6);
%Do the surface fit
options=optimoptions(@lsqcurvefit,'SpecifyObjectiveGradient',false,'MaxFunctionEvaluations',10000);
xdata={r,eta};
params0=linspace(0.1, 0.1, 50);
[params, resnorm]=lsqcurvefit(@modelfun,params0,xdata,g,[],[],options)
for i=1:50
C(i)=params(i);
end
xmin = 1; xmax = 2; dx = 0.01;
etamin=0.10472; etamax=0.55; deta=0.01;
[Xgrid,etagrid]=ndgrid(xmin:dx:xmax,etamin:deta:etamax);
surf(Xgrid,etagrid,modelfun(params,{Xgrid,etagrid}))
zlim([0 8]);
hold on;
scatter3(r(:),eta(:),g(:),'MarkerEdgeColor','none',...
'MarkerFaceColor','b','MarkerFaceAlpha',.5); hold off
xlabel 'x', ylabel '\eta', zlabel 'g(x,\eta)'
function [out,Jacobian]=modelfun(params,xdata)
X=xdata{1};
Y=xdata{2};
for i=1:50
C(i)=params(i);
end
A1 = -0.4194;
A2 = 0.5812;
A3 = 0.6439;
A4 = 0.4730;
eta_c = 0.6452;
PV = 1 + 3*Y./(eta_c-Y)+ A1*(Y./eta_c) + 2*A2*(Y./eta_c).^2 + 3*A3*(Y./eta_c).^3 + 4*A4*(Y./eta_c).^4;
rdf_contact = (PV - 1)./(4*Y);
poly_guess = polyVal2D(C,X,Y,5,5);
out = (poly_guess.*rdf_contact);
if nargout>1
%Jacobian=[X(:), Y(:), ones(size(X(:)))];
end
end

Sign in to comment.


Matt J
Matt J on 20 Mar 2019
Just to sum things up here, below is my implementation of the fit, combining all my various recommendations. Both the algebraic method that I proposed and the method using lsqcurvefit are implemented and compared. You will see that the algebraic method is both faster and more accurate (gives a lower resnorm) than lsqcurvefit.
%% Data Set-up
num=xlsread('example.xlsx',1,'A2:F18110');
Np=4; %polynomial order
g = num(:,6);
otherStuff.r = num(:,4);
otherStuff.eta = num(:,3);
otherStuff.g = g;
otherStuff.eta_c = 0.6452;
otherStuff.Np=Np;
otherStuff.A1 = -0.4194;
otherStuff.A2 = 0.5812;
otherStuff.A3 = 0.6439;
otherStuff.A4 = 0.4730;
%% Fitting (2 methods)
%fit using matrix algebra
tic;
A=func2mat(@(p) modelfun(p,otherStuff), ones(Np+1,Np+1));
Cfit1 =A\g(:);
resnorm1=norm(A*Cfit1-g(:));
toc %Elapsed time is 0.124593 seconds.
%fit iteratively with lsqcurvefit
tic;
options=optimoptions(@lsqcurvefit,'MaxFunctionEvaluations',10000);
Cinitial(1:Np+1,1:Np+1)=0.1;
[Cfit2, resnorm2]=lsqcurvefit(@modelfun,Cinitial,otherStuff,g,[],[],options);
toc %Elapsed time is 1.687990 seconds.
%% Display
figure(1); showFit(Cfit1,otherStuff);
figure(2); showFit(Cfit2,otherStuff);
resnorm1,resnorm2,
function ghat=modelfun(C,otherStuff)
r = otherStuff.r;
eta = otherStuff.eta;
eta_c = otherStuff.eta_c;
Np = otherStuff.Np;
A1 = otherStuff.A1;
A2 = otherStuff.A2;
A3 = otherStuff.A3;
A4 = otherStuff.A4;
PV = 1 + 3*eta./(eta_c-eta)+ A1*(eta./eta_c) + 2*A2*(eta./eta_c).^2 +...
3*A3*(eta./eta_c).^3 + 4*A4*(eta./eta_c).^4;
rdf_contact = (PV - 1)./(4*eta);
Cflip=rot90(reshape(C,Np+1,Np+1),2);
poly_guess = polyVal2D(Cflip,r-1,eta/eta_c,Np,Np);
ghat = (poly_guess.*rdf_contact);
end
function showFit(Cfit,otherStuff)
r = otherStuff.r(:);
eta = otherStuff.eta(:);
g = otherStuff.g(:);
ghat=modelfun(Cfit,otherStuff);
tri = delaunay(r,eta);
%% Plot it with TRISURF
h=trisurf(tri, r, eta, ghat);
h.EdgeColor='b';
h.FaceColor='b';
axis vis3d
hold on;
scatter3(r,eta,g,'MarkerEdgeColor','none',...
'MarkerFaceColor','r','MarkerFaceAlpha',.05);
xlabel 'r', ylabel '\eta', zlabel 'g(r,\eta)'
hold off
end
Both methods give a similar looking surface plot:
untitled.png

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!