Tuesday, November 27, 2012

Exercise 1-7. Write a program to print the value of EOF


Kernighan Ritchie - The C Programming Language (Second Edition)
Esercise 1-7
Write a program to print the value of EOF


/*-------------------------------------------------*/
#include <stdio.h>
int main () {
       
    printf("\nValue of EOF: %d\n", EOF);
   
    return 0;
}
/*-------------------------------------------------*/

The result in my case is -1 .

NOTE
To insert an End Of File (EOF) on a Mac console you need to press "CTRL + D"

Friday, November 9, 2012

How to fit IMPLICIT functions with Matlab.

In Matlab you can perform implicit curve fit using either "nlinfit" or "lsqcurvefit" (the last one only if you have the optimization toolbox) but you need to re-write the implicit function into an explicit one solving it numerically (i.e. using "fzero").
THERE IS NO WAY to fit directly using an implicit expression!

To learn how to successfully fit using an implicit function you can follow the instructions on this Matlab page page or to read and use my worked-out example.

I will show you how my implementation of the Langevin function with a field Weiss (that is an implicit equation). This example is interesting by itself becasue it addresses the topic of how to model an hysteresis loop.
The equation is taken from the paper by D.C. Jiles and D.L. AthertonFerromagnetic Hysteresis, IEEE Transaction on Magnetics, Vol. Mag-19, No. 5, September 1983.

For a general reference about the Langevin equation (and generally for magnetism) I recommend the nice (and inexpensive!) book by S. BlundellMagnetism in Condensed Matter.

So let's start!

So we want to create a fitting function in Matlab using this IMPLICIT equation:




As you see the above equation is implicit due to the term alpha*M.

The first step is to transform the equation into an EXPLICIT one. As I wrote in the beginning, it is necessary just to use the fzero function.
I named this function  LangevinWithWeissField.m (if you are lazy, you can download the file .m directly from this link)

%---------------------------------------------------------------------------------------
function y=LangevinWithWeissField(p, x)% Implicit function: Langevin function with a Weiss field.
% y = A*(coth(B*x(i)+C*y) - 1./(B*x(i)+C*y))
% p is the parameter vector
% Physical meaning of the parameters
% x is the magnetic field H
% y is the magnetization M corresponding to n*<m>z where n is the number of
% magnetic moments,  and <m>z is the projected magnetic moment along the
% field direction (called z) of the total moment mu
%A = amplitude
%B = mu0 m/(kB*T) where kB is the Boltzman constant and T is the temperature
%C is a factor giving the Weiss field (that is proportional to the total
%magnetization y, i.e. M=n*<m>z, note that this term includes also the previous factor
% mu0 m /(kB*T)

% assign the parameters ...

A=p(1);
B=p(2);
C=p(3);

y=zeros(size(x)); % define a vector to allocate the magnetization values
NN=length(x); % total length of the field vector x, i.e. B

opt = optimset('display','off');
% I out off all the messages coming from fzero. If something goes wrong, change this option to 'off'
% to see at which x values fzero failed.

for i=1:NN   
   y(i)=fsolve(@(y) y - A*(coth(B*x(i)+C*y) - 1./(B*x(i)+C*y)) , 0.0001, opt);
   % Here 0.0001 is our starting point to find the solution around 0.
end

end % close the function
%---------------------------------------------------------------------------------------


Once you have saved the code in your Matlab working directory, all the work it's done!

You can plot the Langevin function with the parameter you like, for example (from the Matlab command window)

xx=linspace(-30, 30, 200); % magnetic field
p=[10, 2, 0.1]; %parameters
yy=LangevinWithWeissField(p, xx); % magnetization values
plot(xx, yy, '*-r'); % plot

Now let's see how to "fit" the y-data we just created, i.e. the vector yy.
I use the  LSQCURVEFIT function as done in this mathworks' link.
The code to use in the command line is

lsqcurvefit(@(params, xdata) LangevinWithWeissField(params, xdata),[1 1 1], xx, yy)

In the previous line, the vector  [1 1 1] represents the initial guess of the fitting parameters.
After a while (on my computer it took 54 seconds) the lsqcurvefit gives you the optimal fit parameters, that corresponds to the original values stored in the p vector, namely [10 2 0.1].


To conclude, here is two things to have to keep in mind.
1) It's VERY important to use a good set of values as initial guess for your fit.  With bad (=random) combination of values, your fit could not converge!

2) Using implicit functions can very easily make the fitting SLOW, terribly SLOW!
So, try to avoid their use whenever it is possible.

Any question or comment, as always, is welcome!


UPDATE
You could also be interested in these other posts

(A)  How to hold fitting parameters in nlinfit (MATLAB) 
(B)  edufit: a Matlab data fitting interface
 


I hope that you liked this post and if you did, consider to support my blog! ;) 

Thursday, November 8, 2012

How to hold parameters in nlinfit (Matlab)

If you are using nlinfit to fit your data in Matlab, maybe you are interested (as I was) in finding an easy way to hold some of your parameters (by default there no option to do it).
So how to hold, how to fix some of the parameters?

Looking on the Internet I found this "official" solution link coming directly from the Mathworks support.
It consists in re-writing a new version of nlinfit, called nlinfitsome, that fits only some of parameters. 
The solution is pretty nice and clever but it has a drawback.

The reason I use nlinfit over lsqcurvefit is because it makes relatively easy to get the error bars, however with the proposed solution this advantage is lost so I modified nlinfitsome to address this specific problem.

My function, called mod_nlinfit, has exactly the same arguments of nlinfit and gives you as output the parameters and their error bars. I know that nlinfit has a richer output, but for my needs, having the parameters and error bars is enough.

The code is quite straightforward.

NB
The hold parameters will have no uncertainty.

Credits:
Mathworks support


Matlab code:


function [p, r,  err]=mod_nlinfit(fixed, x, y, fun, beta0, varargin)
% [p, err]=mod_nlinfit(fixed, x, y, fun, beta0, varargin)
% p: vector containing the fitting results

% r: residuals
% err: error bars of the parameters
% fixed is a vector of the same length of beta0 and containing logical elements (i.e. true and
% false). True values mean FIXED parameters.
% EXAMPLE
% If you have three parameters, a vector fixed [true true false]
% means that only the third parameter will change



% all the values in beta0 are copied in p however only the free (NOT fixed)
% ones will change

p=beta0;

% fixed parameters will get zero uncertainity

err=zeros(size(beta0));

% extraction of the free parameters

beta1=beta0(~fixed);

[p1, r1, JJ1] = nlinfit(x,y,@nested_fun,beta1,varargin{:});

% getting the error bars for the free parameters
delta=nlparci(p1, r1, JJ1);
err1=(delta(:, 2)- delta(:, 1))/2;

% update the parameters
p(~fixed)=p1;

%residuals
r=r1;

% ... and the error bars
err(~fixed)=err1;

    % Nested function takes just the parameters to be estimated as inputs
    % It inherits the following from the outer function:
    %   fixed = logical index for fixed elements
    %   beta0 = original guess for the parameters

    function yy=nested_fun(beta, x)
       
        b(fixed)=beta0(fixed);
        b(~fixed)=beta;
       
        yy=fun(b, x);

    end
end




NOTE
This code is used in edufit

Thursday, October 11, 2012

How to check if the input is a number (Matlab)

If you are writing a code where the user is asked to give a number, it can be a good idea to check whether or not the inserted value was really a number instead of something else like a letter or some strange character (%, / _ , etc).

The easiest way to do it is asking for a string input and then attempting (via str2num) to make the conversion string to number. If the conversion fails, then the input wasn't a number!


Here a very simple code that does the job.

------------------------------------------------------------------------------------

function [varargout]=CheckIfNumber(s)
%[varargout]=CheckIfNumber(s)
% s is a string
% varargout is the converted number or row vector


num=str2num(s);

if (isempty(num)==0)

    varargout{1}=num;

end

------------------------------------------------------------------------------------

Wednesday, October 3, 2012

How to get a timestamp (date + time) in Matlab

To complete my  Matlab interface for data fitting, I  needed a Timestamp function, namely a function combining date and time into a single string.

This can be useful when you want to save files having the same base "name" (e.g. foo-x11.dat, foo-x12.dat,  etc) but created at different times.

In Matlab the timestamp is obtained with the command

>>datestr(clock, 0)

which output is (at the time of this post!!!)

10-nov-2012 14:10:20

Unfortunately this string cannot be used into a filename due to the presence of colons.
In addition I would like to have a string starting with "year-month-day" because it is a much more convenient  way to name files created at different times.
 

My version of TimeStamp that produces a string in the format :

year-month-day-hours-h-minutes-m-second-s

as for example

2012-11-10-14h10m20s

Here is the code! Enjoy!
And keep following this blog and my other one.  It's totally worth it!

----------------------------------------------------------------

function [s]=TimeStamp
%[s]=TimeStamp
% time stamp in the format
% year-month-day-hours-h-minutes-m-second-s


% Get current time as date vector

t=clock;

% convert date vector into the string: year-month-day
yyyy_mm_dd=datestr(t, 29);

% convert date vector into the string: day-month-year[blank]hours:minutes:seconds
s=datestr(t, 0);

% get index of the blank
index=find(s==' ');

% remove the string that is before the blank
s=s((index+1):end);

% get the indices of the colons ":"
ii=find(s==':');

%change the first colon into "h" and the second one into "m". Add an "s" at the end of the
%string

h_m_s=strcat(  s(1:(ii(1)-1)), 'h', s( (ii(1)+1): (ii(2)-1)), 'm', s( (ii(2)+1):end), 's');

% concate the two strings
s=strcat(yyyy_mm_dd, '-', h_m_s);


---------------------------------------------------------------- 
Alternatively, you can use the very short code provided by the first commenter :)
In this case:

function s=TimeStamp
s=strcat(datestr(clock,'yyyy-mm-dd-HHMM'),'m',datestr(clock,'ss'),'s');


---------------------------------------------------------------- 
Or the code built on the suggestion provided by the another commenter :)


function s=TimeStamp

tmp = strrep(datestr(clock), ':' , '-' );
s= strrep(tmp, ' ' , '-' );
ii=find(s=='-');
s(ii(4))='h';
s(ii(5))='m';
s=strcat(s, 's');

Tuesday, October 2, 2012

Make Matlab listing ONLY subfolders cointained into a folder

From the Matlab Command Window when you type "dir" you get a list of files and folders contained in your working directory.
The question is: how can you get ONLY the folders?

Here is the solution from stackoverflow.

NOTE
In the following function the output "nameFolds" is a cell variable type!!!

------------------------------------------------------------------------------------------

function nameFolds=ListSubfolders(pathFolder)

d = dir(pathFolder);
isub = [d(:).isdir];
nameFolds = {d(isub).name}';
nameFolds(ismember(nameFolds,{'.','..'})) = [];

 -----------------------------------------------------------------------------------------

Credits:
stackoverflow

Monday, May 3, 2010

Best alternatives to Notepad++ on Mac OS X

On Windows Notepad++ is simply  my favorite programming editor!
Unfortunately, it is not available for Mac OS X :(

Here are some alternatives:

1) Editra: my choice! It's small, neat and free! And ... it's available for Windows too.

2) I haven't tried it yet, but also TextWrangler seems great!