package contribution

View: New views

20 Messages — Rating Filter: Alert me

< Prev | 1 - 2 | Next >

package contribution

by roessli bertrand :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

Dear All,

I have written a collection of octave-files which can be used
to obtain crystal and magnetic structures from single-crystal neutron
diffraction.
The code can be downloaded at

http://spectroscopy.web.psi.ch/tasp/tasp_support.html

I would be happy to contribute a package to octave-forge e.g. in the 'extra' section.
Please tell me if this is a possibility.

I attach the main function and the description file.

Thank You and Best Regards,

Bertrand Roessli

---
http://broessli.redbubble.com/

[mufit.m]

## Copyright (C) 2007,2008,2009 Bertrand Roessli <bertrand.roessli@...> ## ## ## MuFit is free software; you can redistribute it and/or ## modify it under the terms of the GNU General Public ## License as published by the Free Software Foundation; ## either version 2, or (at your option) any later version. ## ## Mufit is distributed in the hope that it will be useful, ## but WITHOUT ANY WARRANTY; without even the implied ## warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR ## PURPOSE. See the GNU General Public License for more ## details. ## ## You should have received a copy of the GNU General Public ## License along with Octave; see the file COPYING. If not, ## write to the Free Software Foundation, Inc., 51 Franklin Street, ## Fifth Floor, Boston, MA 02110-1301, USA. ## ## usage: mufit ## ## Author: Bertrand Roessli <bertrand.roessli@...> ## -*- texinfo -*- ## @deftypefn{Function File} mufit ## @deftypefnx{Function File} mufit ## ## @setfilename mufit.info ## @settitle Mufit Manual: ## @setchapternewpage odd ## @c ## ## @node Top, First Chapter, , (dir) ## @comment node-name, next, previous, up ## ## @menu ## * First Chapter:: ## * Second Chapter:: ## @end menu ## ## @node First Chapter, Second Chapter, Top, Top ## @comment node-name, next, previous, up ## ## @chapter First Chapter ## @cindex Introduction ## ## Mufit calculates or fits neutron nuclear/magnetic intensities or neutron polarisation ## matrices. The program has been tested with Octave-3.2.0. ## ## @itemize @bullet ## @item ## The 'Mufit' function does not take any argument but ## uses 3 input files: 1- 'nuc.inp'; 2- 'crys.inp'; 3- 'mag.inp': ## ## @item ## The nuclear model is given in the file 'nuc.inp'; ## The magnetic model is given in the file 'mag.inp'; ## To fit polarimetry data, the file 'nuc.inp' is required in case ## of nuclear-magnetic interaction terms. ## ## @item ## Lists of HKL/Intensities/Standard deviations are given in ## the files 'nuc.inp' and 'crys.inp'. The polarisation matrices are ## given in the file 'mag.inp'. ## ## @item ## To calculate the neutron polarisation matrices ## it is necessary to give a second vector in the scattering plane ## to define the crystal orientation. ## (see mag.inp for an example) ## ## @item ## The coefficients of the magnetic formfactor are given in the ## function 'formfac.m'. You will probably need to add the coefficients ## for your case, as currently only a few examples are implemented. ## ## @item ## The data can be fitted with the Levenberg-Marquard method or ## the Nelder-Mead algorithm. Please note that the size of the initial ## Simplex depends on the number of parameters. Therefore it is recommended ## that if the Simplex method is used the number of cycles increases with ## increasing number of parameters (~1000 cycles for 20 parameters) ## ## @item ## Simulated Annealing has also been implemented. The algorithm is ## described in the function samin.cc. ## ## @item ## The fit functions have been adapted from the 'optim-1.0.6' package at ## www.octaveforge.org ## ## @item ## The fit options should be given in the files 'nuc.inp' and 'mag.inp'. ## See leasrq.m, nmsmax.m and samin.cc for details. ## Documentation can be obtained by typing "help function_name" at the Octave prompt. ## ## @item ## To install mufit you will need to compile __polmat.cc and samin.cc. ## At the Octave prompt, type mkoctfile ***.cc ## ## @item ## If you use the standalone application, just run the mufit executable. ## The scripts, cpp functions and input files should be in the same directory than ## the executable. Please note that Octave needs to be installed even if ## you prefer to use the executable. ## ## @item ## Please report bugs to bertrand.roessli_at_psi.ch ## ## @end itemize ## @node Second Chapter, First Chapter, Top ## @chapter Second Chapter ## @printindex cp ## @contents ## @seealso{Mufit manual} ## @end deftypefn function exit =mufit %---- main function to calculate/fit neutron polarisation matrices with a few other options clear all; warning("off","all"); exit=1; %---- import everything necessary to calculate and fit polarimetry data, see readmu.m [as,bs,cs,aa,bb,cc,param,dp,flag_k0,flag_nmi,kk0,hhkl_AA,hhkl_rlu,AA1,Pi,Pfo,dPfo,nr,P_1,P_2,moment_type,... fit_type,simplex_options,lmfit_options,sim_ann_par,... lat_cell,ion_cell,moment_cell,propvec_cell_rlu,propvec_cell_AA,spgsym_cell,spgtrans_cell,domains_cell]=readmu; disp(' '); cont=menu('--> Mufit-3.4 (bertrand.roessli@...), available options:', 'calculate the magnetic moment configuration', 'calculate the magnetic interaction vectors/ the structure factors', 'calculate the single-crystal magnetic intensities', 'calculate the polarisation matrices using the Blume''s equations', 'calculate the polarisation matrices/structure factors with domains', 'fit the single-crystal nuclear intensities', 'fit the single-crystal magnetic intensities', 'fit the polarisation matrices', 'draw the magnetic structure', 'exit'); if (cont == 1) %---- calculate the magnetic configuration and print out a list with spin directions %---- useful to check magnetic structure against input file %---- calls: lat2cat.m,moment.m [Mcryst2cart,Mrec2cart]=lat2cart(as,bs,cs,aa,bb,cc); tic [M_cell,spins_cart_cell]=moment(lat_cell,ion_cell,moment_cell,propvec_cell_rlu,propvec_cell_AA,spgsym_cell,spgtrans_cell); toc [fidout]=fopen('mag.str','wt'); fprintf(stdout, 'ION TRANSLATION COORDINATES (RMx,RMy,RMz) (IMx,IMy,IMz) M (m_UB)\n'); %write file to plot magnetic structure with DRAWxtl [fid]=fopen('mag.plot','wt'); color=['Blue'; 'Green'; 'Red'; 'Yellow'; 'Black']; radius=0.1; l=1; %end initialise variables for plot fprintf(fid,'title\n'); fprintf(fid,'cell %f %f %f %f %f %f \n',as,bs,cs,aa,bb,cc); fprintf(fid,'spgp P 1\n'); for k=1:M_cell{1} at_rec_units=Mcryst2cart^-1*M_cell{k+1}.mat'; at_rec_units=[at_rec_units(1)./as,at_rec_units(2)./bs,at_rec_units(3)./cs]; MM_p=Mcryst2cart*[M_cell{k+1}.mpx;M_cell{k+1}.mpy;M_cell{k+1}.mpz]; MM_q=Mcryst2cart*[M_cell{k+1}.mqx;M_cell{k+1}.mqy;M_cell{k+1}.mqz]; MM=sqrt(dot(MM_p,MM_p).+dot(MM_q,MM_q)); %MM=sqrt(norm([M_cell{k}.mpx;M_cell{k}.mpy;M_cell{k}.mpz])^2+norm([M_cell{k}.mqx;M_cell{k}.mqy;M_cell{k}.mqz])^2); fprintf(stdout,'%s: (%+2.2f, %+2.2f, %+2.2f) | (%+2.2f, %+2.2f, %+2.2f)| (%+2.2f, %+2.2f, %+2.2f) | (%+2.2f, %+2.2f, %+2.2f) | %3.2f \n',... M_cell{k+1}.label,M_cell{k+1}.T,at_rec_units,... M_cell{k+1}.mpx,M_cell{k+1}.mpy,M_cell{k+1}.mpz,... M_cell{k+1}.mqx,M_cell{k+1}.mqy,M_cell{k+1}.mqz,... MM); fprintf(fidout,'%s: (%+2.2f, %+2.2f, %+2.2f)| (%+2.2f, %+2.2f, %+2.2f) | (%+2.2f, %+2.2f, %+2.2f) | %3.2f \n',... M_cell{k+1}.label,at_rec_units,... M_cell{k+1}.mpx,M_cell{k+1}.mpy,M_cell{k+1}.mpz,... M_cell{k+1}.mqx,M_cell{k+1}.mqy,M_cell{k+1}.mqz,... MM); fprintf(fid,'atom %s 1 %2.3f %2.3f %2.3f\n',M_cell{k+1}.label,at_rec_units); endfor fclose(fidout); %enter colors for ions flag_sphere=[M_cell{2}.label](1:2); kk=1; while (1) if (kk==1) break; endif; if ([M_cell{kk+1}.label](1:2) ~= flag_sphere) l=l+1; flag_sphere=[M_cell{k+1}.label](1:2); endif if (kk > moment_cell{1}) break;endif; fprintf(fid,'sphere %s %f %s\n',M_cell{kk+1}.label,l*radius,color(l,:)); kk=kk+1; endwhile for k=1:M_cell{1} at_rec_units=Mcryst2cart^-1*M_cell{k+1}.mat'; at_rec_units=[at_rec_units(1)./as,at_rec_units(2)./bs,at_rec_units(3)./cs]; MM_p=Mcryst2cart*[M_cell{k+1}.mpx;M_cell{k+1}.mpy;M_cell{k+1}.mpz]; MM_q=Mcryst2cart*[M_cell{k+1}.mqx;M_cell{k+1}.mqy;M_cell{k+1}.mqz]; MM=sqrt(dot(MM_p,MM_p).+dot(MM_q,MM_q)); sx=(M_cell{k+1}.mpx.+M_cell{k+1}.mqx)/MM; sy=(M_cell{k+1}.mpy.+M_cell{k+1}.mqy)/MM; sz=(M_cell{k+1}.mpz.+M_cell{k+1}.mqz)/MM; if (isnan(sx)) sx=0;endif if (isnan(sy)) sy=0;endif if (isnan(sz)) sz=0;endif fprintf(fid,'arrow %2.3f %2.3f %2.3f %2.3f %2.3f %2.3f %2.3f 0.150 black\n',at_rec_units,sx,sy,sz,MM/2); endfor fclose(fid); fprintf(stdout,'\n-->file mag.str written (input to Drawxtl)\n'); elseif (cont == 2) %---- calculate the magnetic interaction vector and structure factor %---- the magnetic interaction vector is printed out in both crystal and mupad coordinate systems. %---- calls: lat2cart.m,moment.m, magfac.m tic fprintf( stdout,"\nType D H K L Magnetic interaction vector in crystal / MuPAD coordinate systems, F_mag^2:\n"); fprintf(stdout,"-----------------------------------------------------------------------------------------------------------------------------------------------------\n"); [Mcryst2cart,Mrec2cart]=lat2cart(as,bs,cs,aa,bb,cc); [M_cell,spins_cart_cell]=moment(lat_cell,ion_cell,moment_cell,propvec_cell_rlu,propvec_cell_AA,spgsym_cell,spgtrans_cell); ll=1; for kk=1:6:nr for ii=1:domains_cell{1} %---- rotation matrices for domains are given in a cartesian coordinates in input file hkl_all_rlu(1:3,ll)=(Mcryst2cart^-1*domains_cell{ii+1}.rot*Mcryst2cart)*(hhkl_rlu(kk,1:3))'; hkl_all(1:3,ll)=domains_cell{ii+1}.rot*(hhkl_AA(kk,1:3))'; AA1_all(1:3,ll)=domains_cell{ii+1}.rot*AA1(kk,1:3)'; det_r(ll)=domains_cell{ii+1}.det; dtype(ll)=domains_cell{ii+1}.type; Pi_all(1:3,ll)=Pi(kk,1:3)'; popul(ll)=domains_cell{ii+1}.pop; ll++; endfor endfor %---- loop over hkl and calculates polarisation for k=1:ll-1 %Q-vector H=hkl_all(1,k); K=hkl_all(2,k); L=hkl_all(3,k); %----- Q vector in cartesian coordinates (1/A) Q=[H;K;L]'; %---- mag. inter. vec., |F|^2 and magnetic interaction vector in cartesian coordinates [mag_inter_vec,fsqr,mag_vec_cart]=magfac(spins_cart_cell,Q,hkl_all_rlu(:,k)',propvec_cell_rlu,flag_k0,dtype(k)); %---- transform magnetic interaction vector into Cryopad/Mupad coordinate system V1=[H;K;L]; %---- component of vector in scattering plane in cartesian coordinates (1/A) V2=AA1_all(:,k); if (dot(V1/norm(V1),V2/norm(V2)) == 1 ) error("\nredefine scattering plane\n"); return; endif %----- Form unit vectors V1, V2 in scattering plane V3 up V3=cross(V1,V2); V2=cross(V3,V1); V3=V3/sqrt(sum(V3.*V3)); V2=V2/sqrt(sum(V2.*V2)); V1=V1/sqrt(sum(V1.*V1)); V3=V3.*det_r(k); %----- coordinates of HKL in V1, V2, V3 frame %----- transformation matrix U=[V1';V2';V3']; %----- rotate magnetic interaction vector into MuPAD system mag_inter_vec_mu=U*mag_inter_vec; bragg=Mrec2cart^-1*[H;K;L]; %----- print results fprintf(stdout,"%s %+2.3f %+2.3f %+2.3f %+2.3f | (%+2.3f %+2.3f*I, %+2.3f %+2.3f*I, %+2.3f %+2.3f*I) | (%+2.3f %+2.3f*I, %+2.3f %+2.3f*I, %+2.3f %+2.3f*I) | %2.3f\n",... dtype(k),det_r(k),... bragg(1),bragg(2),bragg(3),... real(mag_inter_vec(1)),imag(mag_inter_vec(1)),... real(mag_inter_vec(2)),imag(mag_inter_vec(2)),... real(mag_inter_vec(3)),imag(mag_inter_vec(3)),... real(mag_inter_vec_mu(1)),imag(mag_inter_vec_mu(1)),... real(mag_inter_vec_mu(2)),imag(mag_inter_vec_mu(2)),... real(mag_inter_vec_mu(3)),imag(mag_inter_vec_mu(3)),... fsqr); endfor toc elseif (cont == 3) tic; %-----calculate magnetic single crystal intensities with domains %---- get input from file "crys.inp" (see crys_inp.m); [genpar,singlecrys]=crys_inp(as,bs,cs,aa,bb,cc); %----- get the list of magnetic moments in the unit cell in cartesian coordinates [M_cell,spins_cart_cell]=moment(lat_cell,ion_cell,moment_cell,propvec_cell_rlu,propvec_cell_AA,spgsym_cell,spgtrans_cell); lambda=genpar.lambda; %---- prepare hkl-list with domains, calculate formfactor, Lorenz factor [QQ,fstar,formfactor,popul,llor]=hklprep(as,bs,cs,aa,bb,cc,lambda,singlecrys,domains_cell,spins_cart_cell,propvec_cell_rlu); %----- get magnetic structure factor for the hkl-list with domains [mag_inter_vec,fsqr,mag_vec_cart]=magfac2(spins_cart_cell,fstar,QQ,flag_k0,formfactor); %----- calculate extinction [extinc]=extinc(fsqr',genpar.lambda,genpar.zach,QQ); %----- can proceed with calculations of the intensities %----- magnetic intensity for each hkl if (singlecrys{1}.intfsqr == 0) %---- true if integrated intensities not corrected for Lorentz factor I_hkl=popul.*fsqr'.*llor.*extinc; else I_hkl=popul.*fsqr'.*extinc; endif %----- reshape to sum over the magnetic domains I_hkl_dom=reshape(I_hkl,domains_cell{1},columns(I_hkl)/domains_cell{1})'; %----- calculated magnetic intensities averaged over the domains if (domains_cell{1} > 1) I_crys=sum(I_hkl_dom,2); else I_crys=I_hkl_dom; endif %----- print out calculations nr_hkl=singlecrys{2}; hkl_crys_rlu=[[singlecrys{1,3:nr_hkl.+2}].hhkl_rlu]; fprintf(stdout,'H K L Int Lorenz factor\n'); for k=1:nr_hkl hkl_print=[[singlecrys{1,k+2}].hhkl_rlu]; fprintf(stdout,'%+f %+f %+f %f %f \n',hkl_print(1),hkl_print(2),hkl_print(3),I_crys(k),llor(k*domains_cell{1})); endfor toc elseif (cont == 4) %---- calculate the polarisation matrices according to Blume's formula %---- calls pol.m to calculate scattered polarisation matrix using Blume's equations %---- uses lat2cart.m tic; %---- calculate polarisation matrices with Blume's equations [pf]=pol(reshape(Pi',1,3*nr),hhkl_AA,hhkl_rlu,AA1,nr,flag_k0,flag_nmi,... lat_cell,ion_cell,spgsym_cell,spgtrans_cell,propvec_cell_rlu,propvec_cell_AA,moment_cell); pol_calc=(reshape(pf,3,nr))'; %----- Print out a few information before starting calculations fprintf(stdout,"\nPolarisation matrices calculated from Blume's equations\n"); fprintf(stdout,"! domains are NOT taken into account !\n"); fprintf(stdout,"! efficiency of polariser/analyser assumed to be 1 !\n\n"); fprintf(stdout," H K L || Pol. (incident) | Pol. (measured)| Pol. (calculated)\n"); %----- transform to cartesian system [Mcryst2cart,Mrec2cart]=lat2cart(as,bs,cs,aa,bb,cc); bragg_ind=(Mrec2cart^-1*hhkl_AA')'; fprintf(stdout,'%+2.3f %+2.3f %+2.3f | %+2.2f %+2.2f %+2.2f| %+2.2f %+2.2f %+2.2f | %+2.2f %+2.2f %+2.2f\n',[bragg_ind,Pi,Pfo,pol_calc]'); fprintf(stdout,'\nelapsed time %+3.2f (s)\n',toc); elseif (cont == 5) %---- Calculate polarisation matrices and structure factors^2 with bender efficiencies %---- The polarisation matrices and cross sections are calculated from the neutron cross-section [matrix]=fopen('matrix.lis','wt'); fprintf(stdout,'-->Polarisation matrices and cross-sections for Bender efficiencies of P_1= %+3.2f and P_2= %+3.2f\n\n',P_1,P_2); fprintf(matrix,'-->Polarisation matrices and cross-sections for Bender efficiencies of P_1= %+3.2f and P_2= %+3.2f\n\n',P_1,P_2); %----- transformation matrix to cartesian system for real and reciprocal space [Mcryst2cart,Mrec2cart]=lat2cart(as,bs,cs,aa,bb,cc); [M_cell,spins_cart_cell]=moment(lat_cell,ion_cell,moment_cell,propvec_cell_rlu,propvec_cell_AA,spgsym_cell,spgtrans_cell); ll=1; for kk=1:3:nr for ii=1:domains_cell{1} %---- rotation matrices for domains are given in a cartesian coordinates in input file hkl_all_rlu(1:3,ll)=(Mcryst2cart^-1*domains_cell{ii+1}.rot*Mcryst2cart)*(hhkl_rlu(kk,1:3))'; hkl_all(1:3,ll)=domains_cell{ii+1}.rot*(hhkl_AA(kk,1:3))'; AA1_all(1:3,ll)=domains_cell{ii+1}.rot*AA1(kk,1:3)'; Pi_all(1:3,ll)=Pi(kk,1:3)'; det_r(ll)=domains_cell{ii+1}.det; dtype(ll)=domains_cell{ii+1}.type; popul(ll)=domains_cell{ii+1}.pop; ll++; endfor endfor tic; jj=1; for k=1:domains_cell{1}:ll-1 fprintf(stdout," Final Polarisation Intensity (+ :--> +) Intensity (+ :--> -)\n"); fprintf(stdout,"Type Popul. | Q (rlu) | XX XY XZ | YX YY YZ | ZX ZY ZZ | XX XY XZ | YX YY YZ | ZX ZY ZZ | X-X X-Y X-Z | Y-X Y-Y Y-Z | Z-X Z-Y Z-Z |\n"); fprintf(matrix," Final Polarisation Intensity (+ :--> +) Intensity (+ :--> -)\n"); fprintf(matrix,"Type Popul. | Q (rlu) | XX XY XZ | YX YY YZ | ZX ZY ZZ | XX XY XZ | YX YY YZ | ZX ZY ZZ | X-X X-Y X-Z | Y-X Y-Y Y-Z | Z-X Z-Y Z-Z |\n"); %---- Pi along x I_pp=zeros(3,3); I_pm=zeros(3,3); polar=zeros(3,3); for ii=1:domains_cell{1} H=hkl_all(1,jj); K=hkl_all(2,jj); L=hkl_all(3,jj); %----- Q vector in cartesian coordinates (1/A) Q=[H;K;L]'; Q_rlu=[hkl_all_rlu(1,jj);hkl_all_rlu(2,jj);hkl_all_rlu(3,jj)]'; if (ii == 1) Q_rlu_print=Q_rlu; endif; %---- second vector in scattering plane %---- first vector is given by Q A1=AA1_all(1:3,jj); [I_pp_1d,I_pm_1d,polar_1d]=... feval("onedomain",P_1,P_2,Q,Q_rlu,A1',flag_k0,flag_nmi,spins_cart_cell,dtype(jj),'pi_up',propvec_cell_rlu,det_r(jj)); fprintf(stdout,'%s %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f\n',... dtype(jj),popul(jj),... Q_rlu(1),Q_rlu(2),Q_rlu(3),... polar_1d(1,1),polar_1d(1,2),polar_1d(1,3),polar_1d(2,1),polar_1d(2,2),polar_1d(2,3),polar_1d(3,1),polar_1d(3,2),polar_1d(3,3),... real(I_pp_1d(1,1)), real(I_pp_1d(1,2)),real(I_pp_1d(1,3)),... real(I_pp_1d(2,1)), real(I_pp_1d(2,2)),real(I_pp_1d(2,3)),... real(I_pp_1d(3,1)), real(I_pp_1d(3,2)),real(I_pp_1d(3,3)),... real(I_pm_1d(1,1)), real(I_pm_1d(1,2)),real(I_pm_1d(1,3)),... real(I_pm_1d(2,1)), real(I_pm_1d(2,2)),real(I_pm_1d(2,3)),... real(I_pm_1d(3,1)), real(I_pm_1d(3,2)),real(I_pm_1d(3,3))); fprintf(matrix,'%s %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f\n',... dtype(jj),popul(jj),... Q_rlu(1),Q_rlu(2),Q_rlu(3),... polar_1d(1,1),polar_1d(1,2),polar_1d(1,3),polar_1d(2,1),polar_1d(2,2),polar_1d(2,3),polar_1d(3,1),polar_1d(3,2),polar_1d(3,3),... real(I_pp_1d(1,1)), real(I_pp_1d(1,2)),real(I_pp_1d(1,3)),... real(I_pp_1d(2,1)), real(I_pp_1d(2,2)),real(I_pp_1d(2,3)),... real(I_pp_1d(3,1)), real(I_pp_1d(3,2)),real(I_pp_1d(3,3)),... real(I_pm_1d(1,1)), real(I_pm_1d(1,2)),real(I_pm_1d(1,3)),... real(I_pm_1d(2,1)), real(I_pm_1d(2,2)),real(I_pm_1d(2,3)),... real(I_pm_1d(3,1)), real(I_pm_1d(3,2)),real(I_pm_1d(3,3))); %---- compute intensity and polarisation in each domain I_pp=I_pp.+I_pp_1d*popul(jj); I_pm=I_pm.+I_pm_1d*popul(jj); jj=jj+1; endfor %---- calculate final polarisation polar=(I_pp.-I_pm)./(I_pp.+I_pm); fprintf(stdout,'\naverage for %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f\n\n',... Q_rlu_print(1),Q_rlu_print(2),Q_rlu_print(3),... polar(1,1),polar(1,2),polar(1,3),polar(2,1),polar(2,2),polar(2,3),polar(3,1),polar(3,2),polar(3,3),... real(I_pp(1,1)), real(I_pp(1,2)),real(I_pp(1,3)),... real(I_pp(2,1)), real(I_pp(2,2)),real(I_pp(2,3)),... real(I_pp(3,1)), real(I_pp(3,2)),real(I_pp(3,3)),... real(I_pm(1,1)), real(I_pm(1,2)),real(I_pm(1,3)),... real(I_pm(2,1)), real(I_pm(2,2)),real(I_pm(2,3)),... real(I_pm(3,1)), real(I_pm(3,2)),real(I_pm(3,3))); fprintf(matrix,'\naverage for %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f\n\n',... Q_rlu_print(1),Q_rlu_print(2),Q_rlu_print(3),... polar(1,1),polar(1,2),polar(1,3),polar(2,1),polar(2,2),polar(2,3),polar(3,1),polar(3,2),polar(3,3),... real(I_pp(1,1)), real(I_pp(1,2)),real(I_pp(1,3)),... real(I_pp(2,1)), real(I_pp(2,2)),real(I_pp(2,3)),... real(I_pp(3,1)), real(I_pp(3,2)),real(I_pp(3,3)),... real(I_pm(1,1)), real(I_pm(1,2)),real(I_pm(1,3)),... real(I_pm(2,1)), real(I_pm(2,2)),real(I_pm(2,3)),... real(I_pm(3,1)), real(I_pm(3,2)),real(I_pm(3,3))); endfor jj=1; for k=1:domains_cell{1}:ll-1 fprintf(stdout," Final Polarisation Intensity (- :--> +) Intensity (- :--> -)\n"); fprintf(stdout,"Type Popul. | Q (rlu) | -XX -XY -XZ | -YX -YY -YZ | -ZX -ZY -ZZ | -XX -XY -XZ | -YX -YY -YZ | -ZX -ZY -ZZ | -X-X -X-Y -X-Z | -Y-X -Y-Y -Y-Z | -Z-X -Z-Y -Z-Z |\n"); fprintf(matrix," Final Polarisation Intensity (- :--> +) Intensity (- :--> -)\n"); fprintf(matrix,"Type Popul. | Q (rlu) | -XX -XY -XZ | -YX -YY -YZ | -ZX -ZY -ZZ | -XX -XY -XZ | -YX -YY -YZ | -ZX -ZY -ZZ | -X-X -X-Y -X-Z | -Y-X -Y-Y -Y-Z | -Z-X -Z-Y -Z-Z |\n"); %---- Pi along -x I_mp=zeros(3,3); I_mm=zeros(3,3); polar=zeros(3,3); popul_nz=0; for ii=1:domains_cell{1} H=hkl_all(1,jj); K=hkl_all(2,jj); L=hkl_all(3,jj); %----- Q vector in cartesian coordinates (1/A) Q=[H;K;L]'; Q_rlu=[hkl_all_rlu(1,jj);hkl_all_rlu(2,jj);hkl_all_rlu(3,jj)]'; if (ii == 1) Q_rlu_print=Q_rlu; endif; %---- second vector in scattering plane %---- first vector is given by Q A1=AA1_all(1:3,jj); [I_mp_1d,I_mm_1d,polar_1d]=... feval("onedomain",P_1,P_2,Q,Q_rlu,A1',flag_k0,flag_nmi,spins_cart_cell,dtype(jj),'pi_down',propvec_cell_rlu,det_r(jj)); fprintf(stdout,'%s %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f\n',... dtype(jj),popul(jj),... Q_rlu(1),Q_rlu(2),Q_rlu(3),... polar_1d(1,1),polar_1d(1,2),polar_1d(1,3),polar_1d(2,1),polar_1d(2,2),polar_1d(2,3),polar_1d(3,1),polar_1d(3,2),polar_1d(3,3),... real(I_mp_1d(1,1)), real(I_mp_1d(1,2)),real(I_mp_1d(1,3)),... real(I_mp_1d(2,1)), real(I_mp_1d(2,2)),real(I_mp_1d(2,3)),... real(I_mp_1d(3,1)), real(I_mp_1d(3,2)),real(I_mp_1d(3,3)),... real(I_mm_1d(1,1)), real(I_mm_1d(1,2)),real(I_mm_1d(1,3)),... real(I_mm_1d(2,1)), real(I_mm_1d(2,2)),real(I_mm_1d(2,3)),... real(I_mm_1d(3,1)), real(I_mm_1d(3,2)),real(I_mm_1d(3,3))); fprintf(matrix,'%s %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f\n',... dtype(k),popul(k),... Q_rlu(1),Q_rlu(2),Q_rlu(3),... polar_1d(1,1),polar_1d(1,2),polar_1d(1,3),polar_1d(2,1),polar_1d(2,2),polar_1d(2,3),polar_1d(3,1),polar_1d(3,2),polar_1d(3,3),... real(I_mp_1d(1,1)), real(I_mp_1d(1,2)),real(I_mp_1d(1,3)),... real(I_mp_1d(2,1)), real(I_mp_1d(2,2)),real(I_mp_1d(2,3)),... real(I_mp_1d(3,1)), real(I_mp_1d(3,2)),real(I_mp_1d(3,3)),... real(I_mm_1d(1,1)), real(I_mm_1d(1,2)),real(I_mm_1d(1,3)),... real(I_mm_1d(2,1)), real(I_mm_1d(2,2)),real(I_mm_1d(2,3)),... real(I_mm_1d(3,1)), real(I_mm_1d(3,2)),real(I_mm_1d(3,3))); I_mp=I_mp.+I_mp_1d.*popul(jj); I_mm=I_mm.+I_mm_1d.*popul(jj); jj++; endfor %---- calculate final polarisation polar=(I_mp.-I_mm)./(I_mp.+I_mm); fprintf(stdout,'\naverage for %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f\n\n',... Q_rlu_print(1),Q_rlu_print(2),Q_rlu_print(3),... polar(1,1),polar(1,2),polar(1,3),polar(2,1),polar(2,2),polar(2,3),polar(3,1),polar(3,2),polar(3,3),... real(I_mp(1,1)), real(I_mp(1,2)),real(I_mp(1,3)),... real(I_mp(2,1)), real(I_mp(2,2)),real(I_mp(2,3)),... real(I_mp(3,1)), real(I_mp(3,2)),real(I_mp(3,3)),... real(I_mm(1,1)), real(I_mm(1,2)),real(I_mm(1,3)),... real(I_mm(2,1)), real(I_mm(2,2)),real(I_mm(2,3)),... real(I_mm(3,1)), real(I_mm(3,2)),real(I_mm(3,3))); fprintf(matrix,'\naverage for %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f\n\n',... Q_rlu_print(1),Q_rlu_print(2),Q_rlu_print(3),... polar(1,1),polar(1,2),polar(1,3),polar(2,1),polar(2,2),polar(2,3),polar(3,1),polar(3,2),polar(3,3),... real(I_mp(1,1)), real(I_mp(1,2)),real(I_mp(1,3)),... real(I_mp(2,1)), real(I_mp(2,2)),real(I_mp(2,3)),... real(I_mp(3,1)), real(I_mp(3,2)),real(I_mp(3,3)),... real(I_mm(1,1)), real(I_mm(1,2)),real(I_mm(1,3)),... real(I_mm(2,1)), real(I_mm(2,2)),real(I_mm(2,3)),... real(I_mm(3,1)), real(I_mm(3,2)),real(I_mm(3,3))); endfor fprintf(stdout,'\nelapsed time %+3.2f (s)\n',toc); fprintf(matrix,'\nelapsed time %+3.2f (s)\n',toc); fclose(matrix); elseif (cont == 6) %---- fit nuclear intensities/ structure factors %----- get cell parameters etc... from readnuc [genpar,as,bs,cs,aa,bb,cc,ion_cell,spgsym_cell,spgtrans_cell,singlecrys,nparameters,dnparameters,... nfit_type,simplexnuc_options,lmfitnuc_options,nsim_ann_par]=readnuc; %----- prepare hkl,intensities,errors %--- hkl-list nr_hkl=singlecrys{2}+2; hkl_crys_rlu=[[singlecrys{1,3:nr_hkl}].hhkl_rlu]; hkl_crys_rlu_inp=reshape(hkl_crys_rlu,3,columns(hkl_crys_rlu)/3); hkl_crys_AA=[[singlecrys{1,3:nr_hkl}].hhkl_AA]; hkl_crys_AA_inp=reshape(hkl_crys_AA,3,columns(hkl_crys_AA)/3); %--- intensity list int_crys_inp=[[singlecrys{1,3:nr_hkl}].f2obs]; %--- error list err_crys_inp=1./[[singlecrys{1,3:nr_hkl}].f2sd]; %--- switch intensity <-> structure factor intfsqr=singlecrys{1}.intfsqr; nlat.as=as;nlat.bs=bs;nlat.cs=cs;nlat.aa=aa;nlat.bb=bb;nlat.cc=cc; switch nfit_type %---- fit parameters options case 0 %---- Function to minimize F=@nucfacfit2; %---- numerical derivatives dFdp=@dfdp; %---- additional convergence options minstep=zeros(size(nparameters)); maxstep=Inf*ones(size(nparameters)); options=[minstep, maxstep]; %---- verbose verbose=lmfitnuc_options.verbose; %---- exit options stol=lmfitnuc_options.stol; niter=lmfitnuc_options.niter; %---- print out fitting options disp( '* Levenberg-Marquard algorithm' ); disp( '* fit options:'); fprintf(stdout,'--> Exit if tol. <= %f\n',stol); fprintf(stdout,'--> Max. number of iterations: %i\n', niter); if (verbose == 0) disp('--> fit progress not shown'); else disp('--> fit progess shown'); endif tic [f,p,deltap,pout,pout,ChiSq,R_Bragg,kvg,iter,corp,covp,covr,stdresid,Z,r2]= ... leasqr(int_crys_inp, err_crys_inp, nparameters, dnparameters,F,... stol, niter, verbose, dFdp,options,... hkl_crys_rlu_inp',hkl_crys_AA_inp',... genpar,nlat,intfsqr,... ion_cell,spgsym_cell,spgtrans_cell); toc case 1 %---- Nelder-Mead algorithm stol=simplexnuc_options.stol; niter=simplexnuc_options.niter; fval=simplexnuc_options.fval; sides=simplexnuc_options.sides; verbose=simplexnuc_options.verbose; %---- stop fit options, -1 is to minimise cost function (see %---- nmsmax.m) stopit = [stol niter fval sides verbose -1]; %---- print out fitting options disp( '* Nelder-Mead algorithm' ); disp( '* fit options:'); fprintf(stdout,'--> Exit if simplex size <= %f\n',stol); fprintf(stdout,'--> Max. number of iterations: %i\n', niter); fprintf(stdout,'--> Exit if ChiSq <= %f\n', fval); if (sides == 0) disp('--> uses regular simplex'); else disp('--> uses right-angled simplex'); endif if (verbose == 0) disp('--> fit progress not shown'); else disp('--> fit progess shown'); endif disp('* results saved in file savit'); savit='SAVE SAVIT x fmax nf'; %---- start fit [p, deltap, f, ChiSq, nf, pout, dpout, R_Bragg]=nmsmax("chi",nparameters,dnparameters,stopit,'savit',... "nucfacfit2",int_crys_inp,err_crys_inp,... hkl_crys_rlu_inp',hkl_crys_AA_inp',... genpar,nlat,intfsqr,... ion_cell,spgsym_cell,spgtrans_cell); case 2 %---- fit with simulated annealing %---- control ub=nparameters.*(1.+dnparameters.*sign(nparameters)); lb=nparameters.*(1.-dnparameters.*sign(nparameters)); nt=nsim_ann_par.nt; ns=nsim_ann_par.ns; rt=nsim_ann_par.rt; maxevals=nsim_ann_par.maexevals; neps=nsim_ann_par.neps; functol=nsim_ann_par.functol; paramtol=nsim_ann_par.paramtol; verbosity=nsim_ann_par.verbosity; minarg = 1; # position of parameters in control (counter starts at 0) control = {lb, ub, nt, ns, rt, maxevals, neps, functol, paramtol,verbosity, 1 }; %---- print options disp( '* Simulated-Annealing algorithm' ); disp( '* options:'); fprintf(stdout,'---> # of evaluations before T reduction: %i\n',nt); fprintf(stdout,'---> # of iterations between bounds adjustments: %i\n',ns); fprintf(stdout,'---> T reduction factor: %f\n',rt); fprintf(stdout,'---> maximum iterations: %i\n',maxevals); fprintf(stdout,'---> function tolerance: %f\n',functol); fprintf(stdout,'---> parameters tolerance: %f\n',paramtol); fprintf(stdout,'---> verbosity: %i\n',verbosity); fflush(stdout); unfixed=length(find(dp~=0)); %---- argument of function args={nparameters,"nucfacfit2",int_crys_inp,err_crys_inp,... hkl_crys_rlu_inp',hkl_crys_AA_inp',... genpar,nlat,intfsqr,... ion_cell,spgsym_cell,spgtrans_cell}; %---- start simulated annealing [pout, ChiSq, convergence] = samin("chi", args, control); endswitch if ((nfit_type == 0) ||(nfit_type == 1)) f=f'; %---- print out fit results disp( '* Least Square Estimates of Parameters'); disp(' value sigma(+/-)'); disp([ p deltap ]); disp('* Chi squared'); disp(abs(ChiSq)); disp('* R_Bragg'); disp(R_Bragg); %----- Plot Obs. vs. Calc. Intensities hold on; plot([0:1.1*max(f)],[0:1.1*max(f)],"-r"); plot(int_crys_inp',f',"o","markersize",10); axis([0,1.1*max(f),0,1.1*max(f)],"square"); xlabel("observed");ylabel("calculated");title("fit result"); hold off; %----- Print out results to sreen and to file (crysfit.nuc) %write results to file fprintf('\n\n* Results of fit written in file crysfit.nuc\n'); [fid]=fopen('crysfit.nuc','wt'); fprintf(fid,'\n* R_Bragg = %+f\n', R_Bragg); fprintf(fid,'\n* ChiSq = %+f\n', ChiSq); fprintf(fid,'\n* Calculations :\n'); fprintf(fid,'# H K L Calc Obs StdErr Diff\n'); fprintf(fid,'%+2.3f %+2.3f %+2.3f %+10.3f %+10.3f %+10.3f %+10.3f\n',[hkl_crys_rlu_inp(1:3,:)',f',int_crys_inp',1./err_crys_inp',f'.-int_crys_inp']'); fclose(fid); endif elseif (cont == 7) %---- fit magnetic integrated intensities with Levenberg-Marquard (0), %---- Nelder-Mead (1) or Simulated annealing (2) %---- get crystal parameters [genpar,singlecrys]=crys_inp(as,bs,cs,aa,bb,cc); nr_hkl=singlecrys{2}+2; %---- hkl-list hkl_crys=[[singlecrys{1,3:nr_hkl}].hhkl_rlu]; hkl_meas=reshape(hkl_crys,columns(hkl_crys)/3,3); %---- intensity list int_crys_obs=[[singlecrys{1,3:nr_hkl}].f2obs]; %---- error list err_crys_obs=1./[[singlecrys{1,3:nr_hkl}].f2sd]; %---- intensities intfsqr=singlecrys{1}.intfsqr; %----- get the list of magnetic moments in the unit cell in cartesian coordinates [M_cell,spins_cart_cell]=moment(lat_cell,ion_cell,moment_cell,propvec_cell_rlu,propvec_cell_AA,spgsym_cell,spgtrans_cell); %---- prepare complete hkl-list lambda=genpar.lambda; %---- prepare hkl-list with domains, calculate formfactor, Lorenz factor [QQ,fstar,formfactor,popul,llor]=hklprep(as,bs,cs,aa,bb,cc,lambda,singlecrys,domains_cell,spins_cart_cell,propvec_cell_rlu); %----- switch for fit algorithm switch fit_type %---- Levenberg-Marquard case 0 %---- input parameters pin=param; %---- Function to minimize F=@crys_int_fit2; %---- numerical derivatives dFdp=@dfdp; %---- additional convergence options minstep=zeros(size(pin)); maxstep=Inf*ones(size(pin)); options=[minstep, maxstep]; %---- if verbose 1 output of fit results at the end of leasqr.m global verbose; verbose = 2; %---- exit options stol=lmfit_options.stol; niter=lmfit_options.niter; %---- print out fitting options disp( '* Levenberg-Marquard algorithm' ); disp( '* fit options:'); fprintf(stdout,'--> Exit if tol. <= %f\n',stol); fprintf(stdout,'--> Max. number of iterations: %i\n', niter); if (verbose == 0) disp('--> fit progress not shown'); else disp('--> fit progess shown'); endif tic [f,p,deltap,pout,dpout,ChiSq,R_Bragg,kvg,iter,corp,covp,covr,stdresid,Z,r2]=... leasqr(int_crys_obs,err_crys_obs, pin, dp,F,... stol,niter,verbose, dFdp,options,... hkl_meas,QQ,lat_cell,ion_cell,spgsym_cell,spgtrans_cell,moment_cell,propvec_cell_AA,... fstar,flag_k0,formfactor,genpar,intfsqr,llor,domains_cell); %---- Nelder-Mead case 1 %---- simplex fit parameters options stol=simplex_options.stol; niter=simplex_options.niter; fval=simplex_options.fval; sides=simplex_options.sides; verbose=simplex_options.verbose; stopit = [stol niter fval sides verbose -1]; %---- print out fitting options disp( '* Nelder-Mead algorithm' ); disp( '* fit options:'); fprintf(stdout,'--> Exit if simplex size <= %f\n',stol); fprintf(stdout,'--> Max. number of iterations: %i\n', niter); fprintf(stdout,'--> Exit if ChiSq <= %f\n', fval); if (sides == 0) disp('--> uses regular simplex'); else disp('--> uses right-angled simplex'); endif if (verbose == 0) disp('--> fit progress not shown'); else disp('--> fit progess shown'); endif disp('* results saved in file savit'); savit='SAVE SAVIT x fmax nf'; %---- start fit [p, deltap, f, ChiSq,nf,pout, dpout, R_Bragg]=nmsmax("chi",param,dp,stopit,'savit',... "crys_int_fit2",int_crys_obs,err_crys_obs,... hkl_meas,QQ,lat_cell,ion_cell,spgsym_cell,spgtrans_cell,moment_cell,propvec_cell_AA,... fstar,flag_k0,formfactor,genpar,intfsqr,llor,domains_cell); %---- simulated annealing case 2 %---- fit with simulated annealing %---- control ub=param.*(1.+dp.*sign(param)); lb=param.*(1.-dp.*sign(param)); nt=sim_ann_par.nt; ns=sim_ann_par.ns; rt=sim_ann_par.rt; maxevals=sim_ann_par.maexevals; neps=sim_ann_par.neps; functol=sim_ann_par.functol; paramtol=sim_ann_par.paramtol; verbosity=sim_ann_par.verbosity; minarg = 1; # position of parameters in control (counter starts at 0) control = { lb, ub, nt, ns, rt, maxevals, neps, functol, paramtol,verbosity, 1 }; %---- print options disp( '* Simulated-Annealing algorithm' ); disp( '* options:'); fprintf(stdout,'---> # of evaluations before T reduction: %i\n',nt); fprintf(stdout,'---> # of iterations between bounds adjustments: %i\n',ns); fprintf(stdout,'---> T reduction factor: %f\n',rt); fprintf(stdout,'---> maximum iterations: %i\n',maxevals); fprintf(stdout,'---> function tolerance: %f\n',functol); fprintf(stdout,'---> parameters tolerance: %f\n',paramtol); fprintf(stdout,'---> verbosity: %i\n',verbosity); fflush(stdout); unfixed=length(find(dp~=0)); %---- argument of function args={param,"crys_int_fit2",int_crys_obs,err_crys_obs,... hkl_meas,QQ,lat_cell,ion_cell,spgsym_cell,spgtrans_cell,moment_cell,propvec_cell_AA,... fstar,flag_k0,formfactor,genpar,intfsqr,llor,domains_cell}; %---- start simulated annealing [pout, ChiSq, convergence] = samin("chi", args, control); endswitch if ((fit_type == 0) ||(fit_type == 1)) f=f'; %---- print out fit results disp( '* Least Square Estimates of Parameters'); disp(' value sigma(+/-)'); disp([ p deltap ]); disp('* Chi squared'); disp(abs(ChiSq)); disp('* R_Bragg'); disp(R_Bragg); %---- print out observed/calculate fsqr nr_hkl=singlecrys{2}+2; hkl_crys_rlu=[[singlecrys{1,3:nr_hkl}].hhkl_rlu]; hkl_crys_rlu_inp=reshape(hkl_crys_rlu,3,columns(hkl_crys_rlu)/3); int_crys_inp=[[singlecrys{1,3:nr_hkl}].f2obs]; err_crys_inp=1./[[singlecrys{1,3:nr_hkl}].f2sd]; intfsqr=singlecrys{1}.intfsqr; fprintf(stdout,' H K L Calc Obs StdErr Diff.\n'); fprintf(stdout,'%+2.3f %+2.3f %+2.3f %+10.3f %+10.3f %+10.3f %+10.3f\n',[hkl_crys_rlu_inp(1:3,:)',f',int_crys_inp',1./err_crys_inp',f'.-int_crys_inp']'); [fid]=fopen('crysfit.mag','wt'); fprintf(fid,'* Least Squares Estimates of Parameters\n'); fprintf(fid," %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f\n",pout); fprintf(fid,'* Estimates of Parameters Errors\n'); fprintf(fid," %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f\n",dpout); fprintf(fid,'\n* R_Bragg = %+f\n', R_Bragg); fprintf(fid,'\n* ChiSq = %+f\n', abs(ChiSq)); fprintf(fid,'\n* Calculations :\n'); fprintf(fid,'# H K L Calc Obs StdErr Diff.\n'); fprintf(fid,'%+2.3f %+2.3f %+2.3f %+10.3f %+10.3f %+10.3f %+10.3f\n',[hkl_crys_rlu_inp(1:3,:)',f',int_crys_inp',1./err_crys_inp',f'.-int_crys_inp']'); fclose(fid); fprintf('\n\n* Results of fit written in file crysfit.mag\n'); hold on; plot([0:1.1*max(f)],[0:1.1*max(f)],"-r"); axis([0,1.1*max(f),0,1.1*max(f)],"square"); plot(int_crys_inp',f',"o","markersize",10); xlabel("observed");ylabel("calculated");title("fit result"); hold off; else disp('* Chi squared'); disp(ChiSq); endif elseif (cont == 8) %---- fit of the polarisation matrices switch fit_type case 0 %---- Levenberg-Marquardt algorithm tic %---- augment different lists because of domains ll=1; singlecrys=cell(1,nr/3); singlecrys{2}=nr/3; V_bloc=[]; jj=1; for kk=1:3:nr singlecrys{1,jj+2}.hhkl_AA=hhkl_AA(kk,1:3); singlecrys{1,jj+2}.hhkl_rlu=hhkl_rlu(kk,1:3); for ii=1:domains_cell{1} %---- rotation matrices for domains are given in a cartesian coordinates in input file hkl_all(1:3,ll)=domains_cell{ii+1}.rot*(hhkl_AA(kk,1:3))'; AA1_all(1:3,ll)=domains_cell{ii+1}.rot*AA1(kk,1:3)'; det_r=domains_cell{ii+1}.det; V1=hkl_all(1:3,ll); V2=AA1_all(1:3,ll); V3=cross(V1,V2); V2=cross(V3,V1); if (dot(V1/norm(V1),V2/norm(V2)) == 1 ) error("\nredefine scattering plane\n"); return; endif V1=V1./norm(V1); V2=V2./norm(V2); V3=V3./norm(V3); V3=V3.*det_r; V=[V1';V2';V3']; V_bloc=blkdiag(V_bloc,V); ll++; endfor jj=jj+1; endfor ll=1; Pi_all_cell=cell; for kk=1:3:nr-2 for ii=1:domains_cell{1} Pi_all_cell{ll,1}=Pi(kk:kk+2,1:3); ll++; endfor endfor Pi_all=cell2mat(Pi_all_cell); %----- get the list of magnetic moments in the unit cell in cartesian coordinates [M_cell,spins_cart_cell]=moment(lat_cell,ion_cell,moment_cell,propvec_cell_rlu,propvec_cell_AA,spgsym_cell,spgtrans_cell); %---- prepare hkl-list with domains, calculate formfactor, Lorenz factor lambda=1; %---- used to calculate Lorenz factor; not used anyway for polarisation matrices [QQ,fstar,formfactor,popul,llor]=hklprep(as,bs,cs,aa,bb,cc,lambda,singlecrys,domains_cell,spins_cart_cell,propvec_cell_rlu); %---- input parameters pin=param; %---- Function to minimize F=@polar_fit; %---- numerical derivatives dFdp=@dfdp; %---- additional convergence options minstep=zeros(size(pin)); maxstep=Inf*ones(size(pin)); options=[minstep, maxstep]; %---- if verbose 1 output of fit results at the end of leasqr.m verbose=lmfit_options.verbose; %---- exit options stol=lmfit_options.stol; niter=lmfit_options.niter; %---- print out fitting options disp( '* Levenberg-Marquard algorithm' ); disp( '* fit options:'); fprintf(stdout,'--> Exit if tol. <= %f\n',stol); fprintf(stdout,'--> Max. number of iterations: %i\n', niter); if (verbose == 0) disp('--> fit progress not shown'); else disp('--> fit progess shown'); endif tic err_Pfo=1./reshape(dPfo',1,3*nr); [f,p,deltap,pout,dpout,ChiSq,R_Bragg,kvg,iter,corp,covp,covr,stdresid,Z,r2]=... leasqr(reshape(Pfo',1,columns(Pfo)*rows(Pfo)),err_Pfo, pin, dp,F,... stol,niter,verbose, dFdp,options,... reshape(Pi',3*nr,1),Pi_all,lat_cell,ion_cell,spgsym_cell,spgtrans_cell,propvec_cell_AA,moment_cell,... fstar,QQ,flag_k0,formfactor,... P_1,P_2,V_bloc,flag_nmi,domains_cell); %---- print out fit results disp( '* Least Square Estimates of Parameters'); disp(' value sigma(+/-)'); disp([ p deltap ]); disp('* Chi squared'); disp(abs(ChiSq)); toc %---- fitted polarisation matrices pol_fitted=(reshape(f,3,nr))'; fprintf(stdout,"\n hkl,incident, observed and calculated scattered polarisations\n"); [fid_out]=fopen('mufit.out','wt'); fprintf(fid_out,"\n hkl,incident, observed and calculated scattered polarisations\n"); pol_fitted=(reshape(f,3,nr))'; dpfo=round(100.*dPfo); for k=1:nr fprintf(stdout,"%2.3f %2.3f %2.3f || %+2.2f %+2.2f %+2.2f| %+2.2f (%2i) %+2.2f (%2i) %+2.2f (%2i) | %+2.2f %+2.2f %+2.2f\n",... hhkl_rlu(k,1),hhkl_rlu(k,2),hhkl_rlu(k,3),Pi(k,1),Pi(k,2),Pi(k,3),Pfo(k,1),dpfo(k,1),Pfo(k,2),dpfo(k,2),Pfo(k,3),dpfo(k,3),pol_fitted(k,1),pol_fitted(k,2),pol_fitted(k,3)); fprintf(fid_out,"%2.3f %2.3f %2.3f || %+2.2f %+2.2f %+2.2f| %+2.2f (%2i) %+2.2f (%2i) %+2.2f (%2i) | %+2.2f %+2.2f %+2.2f\n",... hhkl_rlu(k,1),hhkl_rlu(k,2),hhkl_rlu(k,3),Pi(k,1),Pi(k,2),Pi(k,3),Pfo(k,1),dpfo(k,1),Pfo(k,2),dpfo(k,2),Pfo(k,3),dpfo(k,3),pol_fitted(k,1),pol_fitted(k,2),pol_fitted(k,3)); endfor fclose(fid_out); case 1 %----- fit the polarisation matrices with Nelder-Mead method (simplex) %---- augment different lists because of domains ll=1; singlecrys=cell(1,nr/3); singlecrys{2}=nr/3; V_bloc=[]; jj=1; for kk=1:3:nr singlecrys{1,jj+2}.hhkl_AA=hhkl_AA(kk,1:3); singlecrys{1,jj+2}.hhkl_rlu=hhkl_rlu(kk,1:3); for ii=1:domains_cell{1} %---- rotation matrices for domains are given in a cartesian coordinates in input file hkl_all(1:3,ll)=domains_cell{ii+1}.rot*(hhkl_AA(kk,1:3))'; AA1_all(1:3,ll)=domains_cell{ii+1}.rot*AA1(kk,1:3)'; det_r=domains_cell{ii+1}.det; V1=hkl_all(1:3,ll); V2=AA1_all(1:3,ll); V3=cross(V1,V2); V2=cross(V3,V1); if (dot(V1/norm(V1),V2/norm(V2)) == 1 ) error("\nredefine scattering plane\n"); return; endif V1=V1./norm(V1); V2=V2./norm(V2); V3=V3./norm(V3); V3=V3.*det_r; V=[V1';V2';V3']; V_bloc=blkdiag(V_bloc,V); ll++; endfor jj=jj+1; endfor ll=1; Pi_all_cell=cell; for kk=1:3:nr-2 for ii=1:domains_cell{1} Pi_all_cell{ll,1}=Pi(kk:kk+2,1:3); ll++; endfor endfor Pi_all=cell2mat(Pi_all_cell); %----- get the list of magnetic moments in the unit cell in cartesian coordinates [M_cell,spins_cart_cell]=moment(lat_cell,ion_cell,moment_cell,propvec_cell_rlu,propvec_cell_AA,spgsym_cell,spgtrans_cell); %---- prepare hkl-list with domains, calculate formfactor, Lorenz factor lambda=1; %---- used to calculate Lorenz factor; not used anyway for polarisation matrices [QQ,fstar,formfactor,popul,llor]=hklprep(as,bs,cs,aa,bb,cc,lambda,singlecrys,domains_cell,spins_cart_cell,propvec_cell_rlu); %---- fit with nmsmax %---- parameters options stol=simplex_options.stol; niter=simplex_options.niter; fval=simplex_options.fval; sides=simplex_options.sides; verbose=simplex_options.verbose; stopit = [stol niter fval sides verbose -1]; %---- print out fitting options disp( '* Nelder-Mead algorithm' ); disp( '* fit options:'); fprintf(stdout,'--> Exit if simplex size <= %f\n',stol); fprintf(stdout,'--> Max. number of iterations: %i\n', niter); fprintf(stdout,'--> Exit if ChiSq <= %f\n', fval); if (sides == 0) disp('--> uses regular simplex'); else disp('--> uses right-angled simplex'); endif if (verbose == 0) disp('--> fit progress not shown'); else disp('--> fit progess shown'); endif disp('* results saved in file savit'); savit='SAVE SAVIT x fmax nf'; %---- start fit err_Pfo=1./reshape(dPfo',1,3*nr); [p, deltap, f, ChiSq,nf,pout,dpout]=nmsmax("chi",param,dp,stopit,'savit',... "polar_fit",reshape(Pfo',1,columns(Pfo)*rows(Pfo)),err_Pfo,... reshape(Pi',3*nr,1),Pi_all,lat_cell,ion_cell,spgsym_cell,spgtrans_cell,propvec_cell_AA,moment_cell,... fstar,QQ,flag_k0,formfactor,... P_1,P_2,V_bloc,flag_nmi,domains_cell); %---- print out fit results disp( '* Least Square Estimates of Parameters'); disp(' value sigma(+/-)'); disp([ p deltap ]); disp('* Chi squared'); disp(abs(ChiSq)); %---- print out observed and calculated matrices pol_fitted=(reshape(f,3,nr))'; fprintf(stdout,"\n hkl,incident, observed and calculated scattered polarisations\n"); [fid_out]=fopen('mufit.out','wt'); fprintf(fid_out,"\n hkl,incident, observed and calculated scattered polarisations\n"); pol_fitted=(reshape(f,3,nr))'; dpfo=round(100.*dPfo); for k=1:nr fprintf(stdout,"%2.3f %2.3f %2.3f || %+2.2f %+2.2f %+2.2f| %+2.2f (%2i) %+2.2f (%2i) %+2.2f (%2i) | %+2.2f %+2.2f %+2.2f\n",... hhkl_rlu(k,1),hhkl_rlu(k,2),hhkl_rlu(k,3),Pi(k,1),Pi(k,2),Pi(k,3),Pfo(k,1),dpfo(k,1),Pfo(k,2),dpfo(k,2),Pfo(k,3),dpfo(k,3),pol_fitted(k,1),pol_fitted(k,2),pol_fitted(k,3)); fprintf(fid_out,"%2.3f %2.3f %2.3f || %+2.2f %+2.2f %+2.2f| %+2.2f (%2i) %+2.2f (%2i) %+2.2f (%2i) | %+2.2f %+2.2f %+2.2f\n",... hhkl_rlu(k,1),hhkl_rlu(k,2),hhkl_rlu(k,3),Pi(k,1),Pi(k,2),Pi(k,3),Pfo(k,1),dpfo(k,1),Pfo(k,2),dpfo(k,2),Pfo(k,3),dpfo(k,3),pol_fitted(k,1),pol_fitted(k,2),pol_fitted(k,3)); endfor fclose(fid_out); case 2 %---- the following will be taken out when Levenberg-Marquard will %---- use C++ file %---- augment different lists because of domains ll=1; singlecrys=cell(1,nr/3); singlecrys{2}=nr/3; V_bloc=[]; jj=1; for kk=1:3:nr singlecrys{1,jj+2}.hhkl_AA=hhkl_AA(kk,1:3); singlecrys{1,jj+2}.hhkl_rlu=hhkl_rlu(kk,1:3); for ii=1:domains_cell{1} %---- rotation matrices for domains are given in a cartesian coordinates in input file hkl_all(1:3,ll)=domains_cell{ii+1}.rot*(hhkl_AA(kk,1:3))'; AA1_all(1:3,ll)=domains_cell{ii+1}.rot*AA1(kk,1:3)'; det_r=domains_cell{ii+1}.det; V1=hkl_all(1:3,ll); V2=AA1_all(1:3,ll); V3=cross(V1,V2); V2=cross(V3,V1); if (dot(V1/norm(V1),V2/norm(V2)) == 1 ) error("\nredefine scattering plane\n"); return; endif V1=V1./norm(V1); V2=V2./norm(V2); V3=V3./norm(V3); V3=V3.*det_r; V=[V1';V2';V3']; V_bloc=blkdiag(V_bloc,V); ll++; endfor jj=jj+1; endfor ll=1; Pi_all_cell=cell; for kk=1:3:nr-2 for ii=1:domains_cell{1} Pi_all_cell{ll,1}=Pi(kk:kk+2,1:3); ll++; endfor endfor Pi_all=cell2mat(Pi_all_cell); %----- get the list of magnetic moments in the unit cell in cartesian coordinates [M_cell,spins_cart_cell]=moment(lat_cell,ion_cell,moment_cell,propvec_cell_rlu,propvec_cell_AA,spgsym_cell,spgtrans_cell); %---- prepare hkl-list with domains, calculate formfactor, Lorenz factor lambda=1; %---- used to calculate Lorenz factor; not used anyway for polarisation matrices [QQ,fstar,formfactor,popul,llor]=hklprep(as,bs,cs,aa,bb,cc,lambda,singlecrys,domains_cell,spins_cart_cell,propvec_cell_rlu); %---- end of what will be taken out %----- Simulated Annealing %---- control ub=param.*(1.+dp.*sign(param)); lb=param.*(1.-dp.*sign(param)); nt=sim_ann_par.nt; ns=sim_ann_par.ns; rt=sim_ann_par.rt; maxevals=sim_ann_par.maexevals; neps=sim_ann_par.neps; functol=sim_ann_par.functol; paramtol=sim_ann_par.paramtol; verbosity=sim_ann_par.verbosity; minarg = 1; # position of parameters in control (counter starts at 0) control = { lb, ub, nt, ns, rt, maxevals, neps, functol, paramtol,verbosity, 1 }; %---- print options disp( '* Simulated-Annealing algorithm' ); disp( '* options:'); fprintf(stdout,'---> # of evaluations before T reduction: %i\n',nt); fprintf(stdout,'---> # of iterations between bounds adjustments: %i\n',ns); fprintf(stdout,'---> T reduction factor: %f\n',rt); fprintf(stdout,'---> maximum iterations: %i\n',maxevals); fprintf(stdout,'---> function tolerance: %f\n',functol); fprintf(stdout,'---> parameters tolerance: %f\n',paramtol); fprintf(stdout,'---> verbosity: %i\n',verbosity); fflush(stdout); %---- start fit err_Pfo=1./reshape(dPfo',1,3*nr); %---- argument of function args={param,"polar_fit",reshape(Pfo',1,columns(Pfo)*rows(Pfo)),err_Pfo,... reshape(Pi',3*nr,1),Pi_all,lat_cell,ion_cell,spgsym_cell,spgtrans_cell,propvec_cell_AA,moment_cell,... fstar,QQ,flag_k0,formfactor,... P_1,P_2,V_bloc,flag_nmi,domains_cell}; %---- start simulated annealing [pout, ChiSq, convergence] = samin("chi", args, control); endswitch if ((fit_type == 0) ||(fit_type == 1)) %---- write matrices to Latexfile disp('* Observed and calculated matrices witten in ''mufit.tex'''); [fid_tex]=fopen('mufit.tex','wt'); cmnd='\'; fprintf(fid_tex,'%sdocumentclass[10pt,a4paper]{article}\n',cmnd); fprintf(fid_tex,'%susepackage{ifpdf}\n',cmnd); fprintf(fid_tex,'%susepackage[utf8]{inputenc}\n',cmnd); fprintf(fid_tex,'%stitle{Polarisation Matrices}\n',cmnd); fprintf(fid_tex,'%sbegin{document}\n',cmnd); fprintf(fid_tex,'%smaketitle\n',cmnd); fprintf(fid_tex,'%sbegin{center}\n',cmnd); fprintf(fid_tex,'%sbegin{tabular}{ccc|ccc|ccc|ccc}\n',cmnd); ret='\\'; for k=1:nr fprintf(fid_tex,"%2.2f & %2.2f & %2.2f & %1i & %1i & %1i & %2.2f (%i) & %2.2f (%i) & %2.2f (%i) & %2.2f & %2.2f & %2.2f %s\n",... hhkl_rlu(k,1),hhkl_rlu(k,2),hhkl_rlu(k,3),... Pi(k,1),Pi(k,2),Pi(k,3),Pfo(k,1),... dpfo(k,1),Pfo(k,2),dpfo(k,2),... Pfo(k,3),dpfo(k,3),pol_fitted(k,1),pol_fitted(k,2),pol_fitted(k,3),ret); endfor fprintf(fid_tex,'%send{tabular}\n',cmnd); fprintf(fid_tex,'%send{center}\n',cmnd); fprintf(fid_tex,'%send{document}\n',cmnd); fclose(fid_tex); endif elseif (cont == 9) %---- transformation to cartesian matrices [Mcryst2cart,Mrec2cart]=lat2cart(as,bs,cs,aa,bb,cc); %---- magnetic moments components [M_cell,spins_cart_cell]=moment(lat_cell,ion_cell,moment_cell,propvec_cell_rlu,propvec_cell_AA,spgsym_cell,spgtrans_cell); %---- atomic position for k=M_cell{1}:-1:1 at_rec_units(1:3,k)=Mcryst2cart^-1*M_cell{k+1}.mat'; at_rec_units(1:3,k)=[at_rec_units(1,k)./as,at_rec_units(2,k)./bs,at_rec_units(3,k)./cs]; endfor %----- draw magnetic structure lat.as=as; lat.bs=bs; lat.cs=cs; nr_spins=moment_cell{1,1}; polar_cart=moment_cell{1,nr_spins+2}.type; mt_draw_mag(Mcryst2cart,lat, ion_cell, M_cell, polar_cart); elseif (cont == 10) exit=0; return; endif %---- write file with output if (cont == 7 | cont == 8 ) [Mcryst2cart,Mrec2cart]=feval("lat2cart",as,bs,cs,aa,bb,cc); fprintf(stdout, '* Input file written in file ''mag.new''\n'); [file_new]=fopen('mag.new','wt'); fprintf(file_new,'!lattice\n'); fprintf(file_new,'%+2.4f %+2.4f %+2.4f %+3.2f %+3.2f %+3.2f\n',as,bs,cs,aa,bb,cc); fprintf(file_new,'!#atoms\n'); fprintf(file_new,'%i\n',ion_cell{1}); fprintf(file_new,'!AT x y z B occ\n'); for i=1:ion_cell{1} fprintf(file_new,'%s %+1.5f %+1.5f %+1.5f %+1.5f %+1.5f\n',... ion_cell{i+1}.ion,... Mcryst2cart^-1*[ion_cell{i+1}.x./as;ion_cell{i+1}.y./bs;ion_cell{i+1}.z./cs],... ion_cell{i+1}.B,ion_cell{i+1}.occ); endfor fprintf(file_new,'! input switch 0: components along axis RMx,RMy,RMz,Imx,Imy,Imz,phase; 1: polar coordinates Rm,Rhi,Rtheta,Im,Iphi,Iptheta,phase\n'); fprintf(file_new,'%i\n',moment_type); if (moment_type == 0) fprintf(file_new,'! F(Q) Rmx Rmy Rmz Imx Imy Imz phase (units of 2*pi)\n'); fprintf(file_new,'%i\n',moment_cell{1}); for i=1:moment_cell{1} fprintf(file_new,'%s %s',moment_cell{i+1}.ion, moment_cell{i+1}.fq); fprintf(file_new,' %+2.3f %+2.3f %+2.3f',moment_cell{i+1}.mx,moment_cell{i+1}.my, moment_cell{i+1}.mz); fprintf(file_new,' %+2.3f %+2.3f %+2.3f',moment_cell{i+1}.imx, moment_cell{i+1}.imy, moment_cell{i+1}.imz); fprintf(file_new,' %+2.3f\n',moment_cell{i+1}.phase); endfor else fprintf(file_new,'! F(Q) Rm Rphi Rth Im Iphi Ith phase (units of 2*pi)\n'); fprintf(file_new,'%i\n',moment_cell{1}); for i=1:moment_cell{1} fprintf(file_new,'%s %s',moment_cell{i+1}.ion,moment_cell{i+1}.fq); fprintf(file_new,' %+3.3f %+3.3f %+3.3f',moment_cell{i+1}.rm,moment_cell{i+1}.rphi,moment_cell{i+1}.rtheta); fprintf(file_new,' %+3.3f %+3.3f %+3.3f',moment_cell{i+1}.im,moment_cell{i+1}.iphi,moment_cell{i+1}.itheta); fprintf(file_new,' %+2.3f\n',moment_cell{i+1}.phase); endfor endif fprintf(file_new,'!# symmetry operations in space group and symmetry elements\n'); fprintf(file_new,'%i\n',spgsym_cell{1}); [min,max]=size(spgsym_cell); for i=min:max-1 fprintf(file_new,' %s %s %s',spgsym_cell{1,i+1}.x,spgsym_cell{1,i+1}.y,spgsym_cell{1,i+1}.z) fprintf(file_new,' %s %s %s',spgsym_cell{1,i+1}.rmx,spgsym_cell{1,i+1}.rmy,spgsym_cell{1,i+1}.rmz); fprintf(file_new,' %s %s %s',spgsym_cell{1,i+1}.imx,spgsym_cell{1,i+1}.imy,spgsym_cell{1,i+1}.imz); fprintf(file_new,' %2.3f\n',spgsym_cell{1,i+1}.phase); endfor fprintf(file_new,'!# number of translations and translations\n'); fprintf(file_new,'%i\n',spgtrans_cell{1}); for i=1:spgtrans_cell{1} fprintf(file_new,' %i %i %i\n',spgtrans_cell{i+1}.x./as,spgtrans_cell{i+1}.y./bs,spgtrans_cell{i+1}.z./cs); endfor fprintf(file_new,'!# propagation vector\n'); for i=1:propvec_cell_rlu{1} fprintf(file_new,'%s %+1.3f %+1.3f %+1.3f\n',propvec_cell_rlu{i+1}.label,... [propvec_cell_rlu{i+1}.kx;propvec_cell_rlu{i+1}.ky;propvec_cell_rlu{i+1}.kz]); endfor fprintf(file_new,'!#spin domains, populations and parameters\n'); fprintf(file_new,'%i\n',domains_cell{1}); fprintf(file_new,'!population and rotation matrices for domains\n'); for i=1:domains_cell{1} fprintf(file_new,'%1.2f', domains_cell{i+1}.pop); fprintf(file_new,' %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f %+3.2f', domains_cell{i+1}.rot'); fprintf(file_new,' %s\n',domains_cell{i+1}.type); endfor %parameters 1-->10 fprintf(file_new,'!c1 c2 c3 c4 c5 c6 c7 c8 c9 c10\n'); fprintf(file_new,'%+6.2f ',pout(1:10)); fprintf(file_new,'\n'); %parameters shifts 1-->10 fprintf(file_new,'!steps dC1 to DC30\n'); fprintf(file_new,'%+6.2f ',dp(1:10)); fprintf(file_new,'\n'); fprintf(file_new,'!c11 c12 c13 c14 c15 c16 c17 c18 c19 c20\n'); fprintf(file_new,'%+6.2f ',pout(11:20)); fprintf(file_new,'\n'); fprintf(file_new,'!steps dC11 to DC20\n'); fprintf(file_new,'%+6.2f ',dp(11:20)); fprintf(file_new,'\n'); fprintf(file_new,'!c21 c22 c23 c24 c25 c26 c27 c28 c29 c30\n'); fprintf(file_new,'%+6.2f ',pout(21:30)); fprintf(file_new,'\n'); fprintf(file_new,'!steps dC1 to DC30\n'); fprintf(file_new,'%+6.2f ',dp(21:30)); fprintf(file_new,'\n'); fprintf(file_new,'!D1 D2 D3 D4 D5 D6 D7 D8 D9 D10\n'); fprintf(file_new,'%+6.2f ',pout(31:40)); fprintf(file_new,'\n'); fprintf(file_new,'!steps dD1 to DD30\n'); fprintf(file_new,'%+6.2f ',dp(31:40)); fprintf(file_new,'\n'); fprintf(file_new,'! flag for calculation k0 + -k0 (1); flag for pure magnetic reflection (0)\n'); fprintf(file_new,'%i %i\n',flag_k0,flag_nmi); fprintf(file_new,'! bender efficiencies 0.5 < P_1, P_2 < 1\n'); fprintf(file_new,'%1.2f %1.2f\n',P_1,P_2); fclose(file_new); endif endfunction

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

DESCRIPTION (348 bytes) Download Attachment

Re: package contribution

by Søren Hauberg :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

Hi,

fre, 31 07 2009 kl. 22:56 -0700, skrev roessli bertrand:

> I have written a collection of octave-files which can be used
> to obtain crystal and magnetic structures from single-crystal neutron
> diffraction.
> The code can be downloaded at
>
> http://spectroscopy.web.psi.ch/tasp/tasp_support.html
>
> I would be happy to contribute a package to octave-forge e.g. in the
> 'extra' section.
> Please tell me if this is a possibility.

We would love to have your package in Octave-Forge. I had a quick look
at the files you attached, and they look fine, except you link to
'www.octaveforge.org' which doesn't exist. You should link to
'octave.sf.net' instead.

I also had a quick look at the package from the web page you linked to.
This looks like a good package, but it seems like you are including
other packages ('optim' and 'vrml') as part of your package. I can
understand why you want to do this, but if you wish for your package to
be part of Octave-Forge you need to strip these out, and instead make
your package depend on the other packages.

Søren

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by roessli bertrand :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

Dear Soren,

Thanks for your positive answer. I changed the link to octave.sf.net.

Concerning your comment :
"I can
understand why you want to do this, but if you wish for your package to
be part of Octave-Forge you need to strip these out, and instead make
your package depend on the other packages."

I have no problem doing this (and I would prefer...) but I had some issues with the 2 packages.

In the vrml-package:
---------------------------------
1.- vrml_transfo.m:60:% if abs (det (r) - 1) > sqrt (eps), r2 = orthogonalize (r);
I think this function does not exist anymore and was replaced by orth:
vrml_transfo.m:61: if abs (det (r) - 1) > sqrt (eps), r2 = orth (r);

2.- vrml_text.m:49: s = read_options (varargin{:}, "op1",op1,"op0",op0, "default",df);
causes an error in my octave version (3.2). I replaced the line with
s = read_options (varargin, "op1",op1,"op0",op0, "default",df);

In the optim package:
-----------------------------------
Mufit uses nmsmax, leasqr and samin.cc and I modified these functions sligthly
to fit my needs;
-As far as I understand the code, both nmsmax and samin.cc do not allow to keep (some) parameters 'fixed' when fitting and then to release them. For example, the starting temperature in samin.cc diverges if the lower and upper bounds = start parameter values.
-Also in nsmax and leasqr I modified the way the cost functioni is called such as it is possible to keep the parameters within 'bounds'.
-I added a varargin variable to leasqr.m to allow pass additional parameters.

These are small changes, but the program will not run with the functions in
the optim-package as it is.

What should I do??

Thank You and best reagrds,

Bertrand Roessli

---
http://broessli.redbubble.com/

--- On Sat, 8/1/09, Søren Hauberg <soren@...> wrote:

From: Søren Hauberg <soren@...>
Subject: Re: [OctDev] package contribution
To: "roessli bertrand" <broessli@...>
Cc: octave-dev@...
Date: Saturday, August 1, 2009, 5:26 PM

Hi,

fre, 31 07 2009 kl. 22:56 -0700, skrev roessli bertrand:

> I have written a collection of octave-files which can be used
> to obtain crystal and magnetic structures from single-crystal neutron
> diffraction.
> The code can be downloaded at
>
> http://spectroscopy.web.psi.ch/tasp/tasp_support.html
>
> I would be happy to contribute a package to octave-forge e.g. in the
> 'extra' section.
> Please tell me if this is a possibility.

We would love to have your package in Octave-Forge. I had a quick look
at the files you attached, and they look fine, except you link to
'www.octaveforge.org' which doesn't exist. You should link to
'octave.sf.net' instead.

I also had a quick look at the package from the web page you linked to.
This looks like a good package, but it seems like you are including
other packages ('optim' and 'vrml') as part of your package. I can
understand why you want to do this, but if you wish for your package to
be part of Octave-Forge you need to strip these out, and instead make
your package depend on the other packages.

Søren

Re: package contribution

by Søren Hauberg :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

Hi,

man, 03 08 2009 kl. 00:34 -0700, skrev roessli bertrand:

> In the vrml-package:
> ---------------------------------
> 1.- vrml_transfo.m:60:% if abs (det (r) - 1) > sqrt (eps), r2 =
> orthogonalize (r);
> I think this function does not exist anymore and was replaced by orth:
> vrml_transfo.m:61: if abs (det (r) - 1) > sqrt (eps), r2 = orth (r);
>
> 2.- vrml_text.m:49: s = read_options (varargin{:},
> "op1",op1,"op0",op0, "default",df);
> causes an error in my octave version (3.2). I replaced the line with
> s = read_options (varargin, "op1",op1,"op0",op0, "default",df);

These sounds like bugs in the VRML.

Etienne, are these bugs in the package?

Søren

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by Søren Hauberg :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

man, 03 08 2009 kl. 00:34 -0700, skrev roessli bertrand:

> In the optim package:
> -----------------------------------
> Mufit uses nmsmax, leasqr and samin.cc and I modified these functions
> sligthly
> to fit my needs;
> -As far as I understand the code, both nmsmax and samin.cc do not
> allow to keep (some) parameters 'fixed' when fitting and then to
> release them. For example, the starting temperature in samin.cc
> diverges if the lower and upper bounds = start parameter values.
> -Also in nsmax and leasqr I modified the way the cost functioni is
> called such as it is possible to keep the parameters within 'bounds'.
> -I added a varargin variable to leasqr.m to allow pass additional
> parameters.
>
> These are small changes, but the program will not run with the
> functions in
> the optim-package as it is.

The changes you have made are they general purpose?

I would love to hear from some users/developers of the 'optim' package
on this issue.

Søren

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by Olaf Till :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

some users comments (disclaimer: only used leasqr, not nmsmax or
samin):

On Mon, Aug 03, 2009 at 11:25:57AM +0200, Søren Hauberg wrote:

> man, 03 08 2009 kl. 00:34 -0700, skrev roessli bertrand:
> > In the optim package:
> > -----------------------------------
> > Mufit uses nmsmax, leasqr and samin.cc and I modified these functions
> > sligthly
> > to fit my needs;
> > -As far as I understand the code, both nmsmax and samin.cc do not
> > allow to keep (some) parameters 'fixed' when fitting and then to
> > release them. For example, the starting temperature in samin.cc
> > diverges if the lower and upper bounds = start parameter values.

Think it _is_ general purpose. Though one could simply leave these
parameters out, this is not always convenient. Having an option to
keep some parameters 'fixed' is surely desirable. (leasqr has such an
option, but IMO it could be a bit more direct, i.e. an index of
parameters which are to be fixed).

> > -Also in nsmax and leasqr I modified the way the cost functioni is
> > called such as it is possible to keep the parameters within
> > 'bounds'.

I don't know if there are theoretical objections against bounds in
Levenberg-Marquardt optimization, but I also ended up introducing
bounds into leasqr for my use; I would think this is general
purpose. I would favor not to simply reset a parameter to its bound,
but to preserve the direction of the attempted step. Also there is the
question if only bounds or also linear constraints should be
introduced. But I would not know how to implement a possibility for
more than _one_ linear constraint. Maybe there is some theory
available for that ...

> > -I added a varargin variable to leasqr.m to allow pass additional
> > parameters.

This should be unnecessary since one can use anonymous functions for
passing additional parameters to the user function. Also, no
additional arguments to leasqr could then be introduced in the future
without possibly breaking existing code.

> >
> > These are small changes, but the program will not run with the
> > functions in
> > the optim-package as it is.
>
> The changes you have made are they general purpose?
>
> I would love to hear from some users/developers of the 'optim' package
> on this issue.
>
> Søren

For leasqr, it is probably worth making some changes to the original
first (a new last argument for bounds) and so eliminating the need to
duplicate code.

Regards, Olaf

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by Olaf Till :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

On Mon, Aug 03, 2009 at 01:07:23PM +0200, Olaf Till wrote:
> I don't know if there are theoretical objections against bounds in
> Levenberg-Marquardt optimization, but I also ended up introducing
> bounds into leasqr for my use; I would think this is general
> purpose.

Sorry, I remember now that there were some problems. Often, a bound is
desirable because some parameters do not produce useful results (of
the user function) outside a certain range (e.g. a division by zero
might occur). In such a case, the bounds should already be respected
during gradient determination, not only later in the 'step'. This
would mean that if a parameter is already at its bound, only a
_one_-sided gradient can be computed for this parameter, but only on
the 'side' away from the bound. Also, in two-sided bounds, the
possible parameter-change for gradient estimation could become small
(too small for some problems?). Since currently gradient estimation in
leasqr is done by different function (dfdp), which knows nothing on
bounds, implementation is not so easy. One could pass the
bounds-argument also to dfdp, but there might exist user-replacements
for dftp (there is an argument to leasqr just to choose this function)
which do not honour the 'bounds' argument ... so people can not expect
bounds to work without changing their dfdp-replacement.

Olaf

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by roessli bertrand :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

"The changes you have made are they general purpose?"

I hope so, but somebody should check.

Bertrand

---
http://broessli.redbubble.com/

--- On Mon, 8/3/09, Søren Hauberg <soren@...> wrote:

From: Søren Hauberg <soren@...>
Subject: Re: [OctDev] package contribution
To: "roessli bertrand" <broessli@...>
Cc: octave-dev@...
Date: Monday, August 3, 2009, 11:25 AM

man, 03 08 2009 kl. 00:34 -0700, skrev roessli bertrand:

> In the optim package:
> -----------------------------------
> Mufit uses nmsmax, leasqr and samin.cc and I modified these functions
> sligthly
> to fit my needs;
> -As far as I understand the code, both nmsmax and samin.cc do not
> allow to keep (some) parameters 'fixed' when fitting and then to
> release them. For example, the starting temperature in samin.cc
> diverges if the lower and upper bounds = start parameter values.
> -Also in nsmax and leasqr I modified the way the cost functioni is
> called such as it is possible to keep the parameters within 'bounds'.
> -I added a varargin variable to leasqr.m to allow pass additional
> parameters.
>
> These are small changes, but the program will not run with the
> functions in
> the optim-package as it is.

The changes you have made are they general purpose?

I would love to hear from some users/developers of the 'optim' package
on this issue.

Søren

Re: package contribution

by roessli bertrand :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

> > -I added a varargin variable to leasqr.m to allow pass additional
> > parameters.

"This should be unnecessary since one can use anonymous functions for
passing additional parameters to the user function. "

Can you provide an example please?

Thank You,

Bertrand

---
http://broessli.redbubble.com/

--- On Mon, 8/3/09, Olaf Till <olaf.till@...> wrote:

From: Olaf Till <olaf.till@...>
Subject: Re: [OctDev] package contribution
To: octave-dev@...
Date: Monday, August 3, 2009, 1:07 PM

some users comments (disclaimer: only used leasqr, not nmsmax or
samin):

On Mon, Aug 03, 2009 at 11:25:57AM +0200, Søren Hauberg wrote:

> man, 03 08 2009 kl. 00:34 -0700, skrev roessli bertrand:
> > In the optim package:
> > -----------------------------------
> > Mufit uses nmsmax, leasqr and samin.cc and I modified these functions
> > sligthly
> > to fit my needs;
> > -As far as I understand the code, both nmsmax and samin.cc do not
> > allow to keep (some) parameters 'fixed' when fitting and then to
> > release them. For example, the starting temperature in samin.cc
> > diverges if the lower and upper bounds = start parameter values.

Think it _is_ general purpose. Though one could simply leave these
parameters out, this is not always convenient. Having an option to
keep some parameters 'fixed' is surely desirable. (leasqr has such an
option, but IMO it could be a bit more direct, i.e. an index of
parameters which are to be fixed).

> > -Also in nsmax and leasqr I modified the way the cost functioni is
> > called such as it is possible to keep the parameters within
> > 'bounds'.

I don't know if there are theoretical objections against bounds in
Levenberg-Marquardt optimization, but I also ended up introducing
bounds into leasqr for my use; I would think this is general
purpose. I would favor not to simply reset a parameter to its bound,
but to preserve the direction of the attempted step. Also there is the
question if only bounds or also linear constraints should be
introduced. But I would not know how to implement a possibility for
more than _one_ linear constraint. Maybe there is some theory
available for that ...

> > -I added a varargin variable to leasqr.m to allow pass additional
> > parameters.

This should be unnecessary since one can use anonymous functions for
passing additional parameters to the user function. Also, no
additional arguments to leasqr could then be introduced in the future
without possibly breaking existing code.

> >
> > These are small changes, but the program will not run with the
> > functions in
> > the optim-package as it is.
>
> The changes you have made are they general purpose?
>
> I would love to hear from some users/developers of the 'optim' package
> on this issue.
>
> Søren

For leasqr, it is probably worth making some changes to the original
first (a new last argument for bounds) and so eliminating the need to
duplicate code.

Regards, Olaf

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by Olaf Till :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

On Mon, Aug 03, 2009 at 06:33:09AM -0700, roessli bertrand wrote:

> > > -I added a varargin variable to leasqr.m to allow pass additional
> > > parameters.
>
> "This should be unnecessary since one can use anonymous functions for
> passing additional parameters to the user function. "
>
> Can you provide an example please?
>
> Thank You,
>
> Bertrand

Here is an example:

octave:1> function ret = f (x, p, additional_p)
> ret = log (x + p(1)) + additional_p.base;
> endfunction
octave:2> [y, p, cvg] = leasqr ((1:10)(:), [1.7; 1.8; 1.9; 1.9; 1.95; 2.0; 2.05; 2.1; 2.05; 2.2], [0], @ (x_anon, p_anon) feval ("f", x_anon, p_anon, struct ("base", 1)))
y =

0.91618
1.65211
2.07144
2.36599
2.59323
2.77827
2.93436
3.06934
3.18825
3.29451

p = -0.080407
cvg = 1
octave:3>

Regards, Olaf

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by Etienne Grossmann-2 :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

Hi Soren, hi Roessli,

I just checked and yes, they are bugs. I made the changes and will
test and commit them tonight.

Thanks,

Etienne

On Mon, Aug 3, 2009 at 2:24 AM, Søren Hauberg<soren@...> wrote:

> Hi,
>
> man, 03 08 2009 kl. 00:34 -0700, skrev roessli bertrand:
>> In the vrml-package:
>> ---------------------------------
>> 1.- vrml_transfo.m:60:% if abs (det (r) - 1) > sqrt (eps), r2 =
>> orthogonalize (r);
>> I think this function does not exist anymore and was replaced by orth:
>> vrml_transfo.m:61: if abs (det (r) - 1) > sqrt (eps), r2 = orth (r);
>>
>> 2.- vrml_text.m:49: s = read_options (varargin{:},
>> "op1",op1,"op0",op0, "default",df);
>> causes an error in my octave version (3.2). I replaced the line with
>> s = read_options (varargin, "op1",op1,"op0",op0, "default",df);
>
> These sounds like bugs in the VRML.
>
> Etienne, are these bugs in the package?
>
> Søren
>
>

Re: package contribution

by Etienne Grossmann-2 :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

Hi All,

thanks Soren, thanks Roessli, I just commited the changes to the bugs
that you pointed out to me.

Cheers,

Etienne

On Mon, Aug 3, 2009 at 8:03 AM, Etienne Grossmann<etienne@...> wrote:

> Hi Soren, hi Roessli,
>
> I just checked and yes, they are bugs. I made the changes and will
> test and commit them tonight.
>
> Thanks,
>
> Etienne
>
> On Mon, Aug 3, 2009 at 2:24 AM, Søren Hauberg<soren@...> wrote:
>> Hi,
>>
>> man, 03 08 2009 kl. 00:34 -0700, skrev roessli bertrand:
>>> In the vrml-package:
>>> ---------------------------------
>>> 1.- vrml_transfo.m:60:% if abs (det (r) - 1) > sqrt (eps), r2 =
>>> orthogonalize (r);
>>> I think this function does not exist anymore and was replaced by orth:
>>> vrml_transfo.m:61: if abs (det (r) - 1) > sqrt (eps), r2 = orth (r);
>>>
>>> 2.- vrml_text.m:49: s = read_options (varargin{:},
>>> "op1",op1,"op0",op0, "default",df);
>>> causes an error in my octave version (3.2). I replaced the line with
>>> s = read_options (varargin, "op1",op1,"op0",op0, "default",df);
>>
>> These sounds like bugs in the VRML.
>>
>> Etienne, are these bugs in the package?
>>
>> Søren
>>
>>
>

Re: package contribution

by roessli bertrand :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

I learned something, Thank You
Bertrand
---
http://broessli.redbubble.com/

--- On Mon, 8/3/09, Olaf Till <olaf.till@...> wrote:

From: Olaf Till <olaf.till@...>
Subject: Re: [OctDev] package contribution
To: octave-dev@...
Date: Monday, August 3, 2009, 4:36 PM

On Mon, Aug 03, 2009 at 06:33:09AM -0700, roessli bertrand wrote:

> > > -I added a varargin variable to leasqr.m to allow pass additional
> > > parameters.
>
> "This should be unnecessary since one can use anonymous functions for
> passing additional parameters to the user function. "
>
> Can you provide an example please?
>
> Thank You,
>
> Bertrand

Here is an example:

octave:1> function ret = f (x, p, additional_p)
> ret = log (x + p(1)) + additional_p.base;
> endfunction
octave:2> [y, p, cvg] = leasqr ((1:10)(:), [1.7; 1.8; 1.9; 1.9; 1.95; 2.0; 2.05; 2.1; 2.05; 2.2], [0], @ (x_anon, p_anon) feval ("f", x_anon, p_anon, struct ("base", 1)))
y =

   0.91618
   1.65211
   2.07144
   2.36599
   2.59323
   2.77827
   2.93436
   3.06934
   3.18825
   3.29451

p = -0.080407
cvg = 1
octave:3>

Regards, Olaf

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by Olaf Till :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

On Mon, Aug 03, 2009 at 01:49:45PM +0200, Olaf Till wrote:

> On Mon, Aug 03, 2009 at 01:07:23PM +0200, Olaf Till wrote:
> > I don't know if there are theoretical objections against bounds in
> > Levenberg-Marquardt optimization, but I also ended up introducing
> > bounds into leasqr for my use; I would think this is general
> > purpose.
>
> Sorry, I remember now that there were some problems. Often, a bound is
> desirable because some parameters do not produce useful results (of
> the user function) outside a certain range (e.g. a division by zero
> might occur). In such a case, the bounds should already be respected
> during gradient determination, not only later in the 'step'. This
> would mean that if a parameter is already at its bound, only a
> _one_-sided gradient can be computed for this parameter, but only on
> the 'side' away from the bound. Also, in two-sided bounds, the
> possible parameter-change for gradient estimation could become small
> (too small for some problems?). Since currently gradient estimation in
> leasqr is done by different function (dfdp), which knows nothing on
> bounds, implementation is not so easy. One could pass the
> bounds-argument also to dfdp, but there might exist user-replacements
> for dftp (there is an argument to leasqr just to choose this function)
> which do not honour the 'bounds' argument ... so people can not expect
> bounds to work without changing their dfdp-replacement.

Some new suggestions.

-- "bounds" in leasqr

Attached are new versions of leasqr.m and dfdp.m. They handle "bounds"
in steps and in gradient estimation. The "options" argument is used
for this and is now a structure of options. For backwards
compatibility, the old options-matrix is still recognized. A lot of
cosmetic changes have been made to make leasqr.m better editable under
emacs, so almost each line is changed and its not much use making
diffs :-(. I have tested it a bit. Since leasqr does not seem to have
an "assigned" active maintainer, I would commit it and change the
version number of the optim package, but Soeren please decide on that.

Roessli, if there is nothing wrong with the new leasqr version I
suggest you rewrite your code to make use of it, and make it depend on
the new version of the optim package. The computations according to
Bragg you had put inside it can surely be made _after_ calling leasqr?

-- "bounds" in nmsmax

I am not familiar with the method in nmsmax, but at first glance I do
not see where you introduced "bounds" ...

-- keeping parameters "fixed"

The way you have done this in nmsmax.m could break existing code since
the new argument is not the last. And it can not be made the last due
to the already existing varargin argument. (BTW, why not using an
argument which is a logical vector, true for each parameter to be
fixed?) To be on the safe side, I would suggest you to solve the
"fixing" problem with a wrapper, which could be left in your package
and need not necessarily put into the optim package. Something like
the following sketch (may have overlooked some details):

function [outp1, outp2, ...] = optim_wrapper (fixed, varargin)
## according to index "fixed", remove fixed parameters and possibly
## other related vector-elements
...
## replace name of (or handle to) user-function with an anonymous
## function like that (calls the here defined subfunction)
@ (x, p) local_func (x, p, user_func, fixed, values_of_fixed_parameters)
...
## call optimization routine
[oup1, outp2, ...] = call_optimizer (varargin(:));
## reintroduce fixed parameters into output arguments if necessary
...
## define subfunction
function ret = local_func (x, p, user_func, fixed, values_fixed)
## reintroduce fixed parameters into "p" according to "fixed" and
## "values_fixed"
...
## call user function
ret = feval (user_func, x, p);

% Copyright (C) 1992-1994 Richard Shrager
% Copyright (C) 1992-1994 Arthur Jutan
% Copyright (C) 1992-1994 Ray Muzic
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; If not, see <http://www.gnu.org/licenses/>.

function prt=dfdp(x,f,p,dp,func,bounds)
% numerical partial derivatives (Jacobian) df/dp for use with leasqr
% --------INPUT VARIABLES---------
% x=vec or matrix of indep var(used as arg to func) x=[x0 x1 ....]
% f=func(x,p) vector initialsed by user before each call to dfdp
% p= vec of current parameter values
% dp= fractional increment of p for numerical derivatives
% dp(j)>0 central differences calculated
% dp(j)<0 one sided differences calculated
% dp(j)=0 sets corresponding partials to zero; i.e. holds p(j) fixed
% func=string naming the function (.m) file
% e.g. to calc Jacobian for function expsum prt=dfdp(x,f,p,dp,'expsum')
% bounds=two-column-matrix of lower and upper bounds for parameters
%----------OUTPUT VARIABLES-------
% prt= Jacobian Matrix prt(i,j)=df(i)/dp(j)
%================================

m=size(x,1); if (m==1), m=size(x,2); endif %# PAK: in case #cols > #rows
n=length(p); %dimensions
prt=zeros(m,n);del=zeros(n,1); % initialise Jacobian to Zero
del = dp .* p; %cal delx=fract(dp)*param value(p)
idx = p == 0;
del(idx) = dp(idx); %if param=0 delx=fraction
idx = dp > 0;
del(idx) = abs (del(idx)); % not for one-sided intervals, changed
% direction of intervals could change
% behavior of optimization without bounds
min_del = min (abs (del), bounds(:, 2) - bounds(:, 1));
for j=1:n
ps = p;
if (dp(j)~=0)
if (dp(j) < 0)
ps(j) = p(j) + del(j);
if (ps(j) < bounds(j, 1) || ps(j) > bounds(j, 2))
t_del1 = max (bounds(j, 1) - p(j), - abs (del(j))); %
%non-positive
t_del2 = min (bounds(j, 2) - p(j), abs (del(j))); %
%non-negative
if (- t_del1 > t_del2)
del(j) = t_del1;
else
del(j) = t_del2;
endif
ps(j) = p(j) + del(j);
endif
prt(:, j) = (feval (func, x, ps) - f) / del(j);
else
if (p(j) - del(j) < bounds(j, 1))
tp = bounds(j, 1);
ps(j) = tp + min_del(j);
elseif (p(j) + del(j) > bounds(j, 2))
ps(j) = bounds(j, 2);
tp = ps(j) - min_del(j);
else
ps(j) = p(j) + del(j);
tp = p(j) - del(j);
min_del(j) = 2 * del(j);
endif
f1 = feval (func, x, ps);
ps(j) = tp;
prt(:, j) = (f1 - feval (func, x, ps)) / (min_del(j));
endif
endif
endfor

% Copyright (C) 1992-1994 Richard Shrager
% Copyright (C) 1992-1994 Arthur Jutan
% Copyright (C) 1992-1994 Ray Muzic
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; If not, see <http://www.gnu.org/licenses/>.

function [f, p, cvg, iter, corp, covp, covr, stdresid, Z, r2] = ...
leasqr (x, y, pin, F, stol, niter, wt, dp, dFdp, options)

%% function [f, p, cvg, iter, corp, covp, covr, stdresid, Z, r2] =
%% leasqr (x, y, pin, F, {stol, niter, wt, dp, dFdp, options})
%%
%% Levenberg-Marquardt nonlinear regression of f(x,p) to y(x).
%%
%% Version 3.beta
%% Optional parameters are in braces {}.
%% x = column vector or matrix of independent variables, 1 row per
%% observation: x = [x0 x1....xm].
%% y = column vector of observed values, same number of rows as x.
%% wt = column vector (dim=length(x)) of statistical weights. These
%% should be set to be proportional to (sqrt of var(y))^-1; (That is,
%% the covariance matrix of the data is assumed to be proportional to
%% diagonal with diagonal equal to (wt.^2)^-1. The constant of
%% proportionality will be estimated.); default = ones(length(y),1).
%% pin = column vec of initial parameters to be adjusted by leasqr.
%% dp = fractional increment of p for numerical partial derivatives;
%% default = .001*ones(size(pin))
%% dp(j) > 0 means central differences on j-th parameter p(j).
%% dp(j) < 0 means one-sided differences on j-th parameter p(j).
%% dp(j) = 0 holds p(j) fixed i.e. leasqr wont change initial guess: pin(j)
%% F = name of function in quotes; the function shall be of the form y=f(x,p),
%% with y, x, p of the form y, x, pin as described above.
%% dFdp = name of partial derivative function in quotes; default is "dfdp", a
%% slow but general partial derivatives function; the function shall be
%% of the form prt=dfdp(x,f,p,dp,F) (see dfdp.m).
%% stol = scalar tolerance on fractional improvement in scalar sum of
%% squares = sum((wt.*(y-f))^2); default stol = .0001;
%% niter = scalar maximum number of iterations; default = 20;
%%
%% options = structure, currently recognized fields are 'fract_prec',
%% 'max_fract_change', and 'bounds'. For backwards compatibility,
%% 'options' can also be a matrix whose first and second column
%% contains the values of 'fract_prec' and 'max_fract_change',
%% respectively.
%%
%% Field 'options.fract_prec': column vector (same length as 'pin') of
%% desired fractional precisions in parameter estimates. Iterations
%% are terminated if change in parameter vector (chg) on two
%% consecutive iterations is less than their corresponding elements in
%% 'options.fract_prec'. [ie. all(abs(chg*current parm est) <
%% options.fract_prec) on two consecutive iterations.], default =
%% zeros().
%%
%% Field 'options.max_fract_change': column vector (same length as
%% 'pin) of maximum fractional step changes in parameter vector.
%% Fractional change in elements of parameter vector is constrained to
%% be at most 'options.max_fract_change' between sucessive iterations.
%% [ie. abs(chg(i))=abs(min([chg(i)
%% options.max_fract_change(i)*current param estimate])).], default =
%% Inf*ones().
%%
%% Field 'options.bounds': two-column-matrix, one row for each
%% parameter in 'pin'. Each row contains a minimal and maximal value
%% for each parameter. Default: [-Inf, Inf] in each row.
%%
%%
%% OUTPUT VARIABLES
%% f = column vector of values computed: f = F(x,p).
%% p = column vector trial or final parameters. i.e, the solution.
%% cvg = scalar: = 1 if convergence, = 0 otherwise.
%% iter = scalar number of iterations used.
%% corp = correlation matrix for parameters.
%% covp = covariance matrix of the parameters.
%% covr = diag(covariance matrix of the residuals).
%% stdresid = standardized residuals.
%% Z = matrix that defines confidence region (see comments in the source).
%% r2 = coefficient of multiple determination.
%%
%% All Zero guesses not acceptable

%% A modified version of Levenberg-Marquardt
%% Non-Linear Regression program previously submitted by R.Schrager.
%% This version corrects an error in that version and also provides
%% an easier to use version with automatic numerical calculation of
%% the Jacobian Matrix. In addition, this version calculates statistics
%% such as correlation, etc....
%%
%% Version 3 Notes
%% Errors in the original version submitted by Shrager (now called version 1)
%% and the improved version of Jutan (now called version 2) have been corrected.
%% Additional features, statistical tests, and documentation have also been
%% included along with an example of usage. BEWARE: Some the the input and
%% output arguments were changed from the previous version.
%%
%% Ray Muzic <rfm2@...>
%% Arthur Jutan <jutan@...>

%% Richard I. Shrager (301)-496-1122
%% Modified by A.Jutan (519)-679-2111
%% Modified by Ray Muzic 14-Jul-1992
%% 1) add maxstep feature for limiting changes in parameter estimates
%% at each step.
%% 2) remove forced columnization of x (x=x(:)) at beginning. x could be
%% a matrix with the ith row of containing values of the
%% independent variables at the ith observation.
%% 3) add verbose option
%% 4) add optional return arguments covp, stdresid, chi2
%% 5) revise estimates of corp, stdev
%% Modified by Ray Muzic 11-Oct-1992
%% 1) revise estimate of Vy. remove chi2, add Z as return values
%% Modified by Ray Muzic 7-Jan-1994
%% 1) Replace ones(x) with a construct that is compatible with versions
%% newer and older than v 4.1.
%% 2) Added global declaration of verbose (needed for newer than v4.x)
%% 3) Replace return value var, the variance of the residuals with covr,
%% the covariance matrix of the residuals.
%% 4) Introduce options as 10th input argument. Include
%% convergence criteria and maxstep in it.
%% 5) Correct calculation of xtx which affects coveraince estimate.
%% 6) Eliminate stdev (estimate of standard deviation of parameter
%% estimates) from the return values. The covp is a much more
%% meaningful expression of precision because it specifies a confidence
%% region in contrast to a confidence interval.. If needed, however,
%% stdev may be calculated as stdev=sqrt(diag(covp)).
%% 7) Change the order of the return values to a more logical order.
%% 8) Change to more efficent algorithm of Bard for selecting epsL.
%% 9) Tighten up memory usage by making use of sparse matrices (if
%% MATLAB version >= 4.0) in computation of covp, corp, stdresid.
%% Modified by Francesco Potortì
%% for use in Octave
%% Added bounds feature to this file and to dfdp.m, Olaf Till 4-Aug-2009
%%
%% References:
%% Bard, Nonlinear Parameter Estimation, Academic Press, 1974.
%% Draper and Smith, Applied Regression Analysis, John Wiley and Sons, 1981.
%%
%%set default args

%% argument processing
%%

%% if (sscanf(version,'%f') >= 4),
vernum = sscanf(version,"%f");
if vernum(1) >= 4,
global verbose
plotcmd="plot(x(:,1),y,''+'',x(:,1),f); figure(gcf)";
else
plotcmd="plot(x(:,1),y,''+'',x(:,1),f); shg";
endif
if (exist("OCTAVE_VERSION"))
global verbose
plotcmd='plot(x(:,1),y,"+;data;",x(:,1),f,";fit;");';
endif

if(exist("verbose")~=1), %If verbose undefined, print nothing
verbose=0; %This will not tell them the results
endif

if (nargin <= 8), dFdp='dfdp'; endif
if (nargin <= 7), dp=.001*(pin*0+1); endif %DT
if (nargin <= 6), wt=ones(length(y),1); endif % SMB modification
if (nargin <= 5), niter=20; endif
if (nargin == 4), stol=.0001; endif
%%

y=y(:); wt=wt(:); pin=pin(:); dp=dp(:); %change all vectors to columns
%% check data vectors- same length?
m=length(y); n=length(pin); [m1,m2]=size(x);
if m1~=m ,error('input(x)/output(y) data must have same number of rows ') ,endif

%% processing of 'options'
pprec = zeros (n, 1);
maxstep = Inf * ones (n, 1);
bounds = cat (2, -Inf * ones (n, 1), Inf * ones (n, 1));
if (nargin > 9)
if (ismatrix (options)) % backwards compatibility
tp = options;
options = struct ("fract_prec", tp(:, 1));
if (columns (tp) > 1)
options.max_fract_change = tp(:, 2);
endif
endif
if (isfield (options, "fract_prec"))
pprec = options.fract_prec;
if (rows (pprec) != n || columns (pprec) != 1)
error ("fractional precisions: wrong dimensions");
endif
endif
if (isfield (options, "max_fract_change"))
maxstep = options.max_fract_change;
if (rows (maxstep) != n || columns (maxstep) != 1)
error ("maximum fractional step changes: wrong dimensions");
endif
endif
if (isfield (options, "bounds"))
bounds = options.bounds;
if (rows (bounds) != n || columns (bounds) != 2)
error ("bounds: wrong dimensions");
endif
idx = bounds(:, 1) > bounds(:, 2);
tp = bounds(idx, 2);
bounds(idx, 2) = bounds(idx, 1);
bounds(idx, 1) = tp;
if (any (idx = bounds(:, 1) == bounds(:, 2)))
warning ("lower and upper bounds identical for some parameters, setting the respective elements of 'dp' to zero");
dp(idx) = 0;
endif
if (any (idx = pin < bounds(:, 1)))
warning ("some initial parameters set to lower bound");
pin(idx) = bounds(idx, 1);
endif
if (any (idx = pin > bounds(:, 2)))
warning ("some initial parameters set to upper bound");
pin(idx) = bounds(idx, 2);
endif
endif
endif

if (all (dp == 0))
error ("no free parameters");
endif

%% set up for iterations
p = pin;
f=feval(F,x,p); fbest=f; pbest=p;
r=wt.*(y-f);
ss=r'*r;
sbest=ss;
nrm=zeros(n,1);
chgprev=Inf*ones(n,1);
cvg=0;
epsLlast=1;
epstab=[.1, 1, 1e2, 1e4, 1e6];

%% do iterations
for iter=1:niter,
pprev=pbest;
prt=feval(dFdp,x,fbest,pprev,dp,F,bounds);
r=wt.*(y-fbest);
sprev=sbest;
sgoal=(1-stol)*sprev;
for j=1:n,
if dp(j)==0,
nrm(j)=0;
else
prt(:,j)=wt.*prt(:,j);
nrm(j)=prt(:,j)'*prt(:,j);
if nrm(j)>0,
nrm(j)=1/sqrt(nrm(j));
endif
endif
prt(:,j)=nrm(j)*prt(:,j);
endfor
%% above loop could ? be replaced by:
%% prt=prt.*wt(:,ones(1,n));
%% nrm=dp./sqrt(diag(prt'*prt));
%% prt=prt.*nrm(:,ones(1,m))';
[prt,s,v]=svd(prt,0);
s=diag(s);
g=prt'*r;
for jjj=1:length(epstab),
epsL = max(epsLlast*epstab(jjj),1e-7);
se=sqrt((s.*s)+epsL);
gse=g./se;
chg=((v*gse).*nrm);
%% check the change constraints and apply as necessary
ochg=chg;
idx = ~isinf(maxstep);
limit = abs(maxstep(idx).*pprev(idx));
chg(idx) = min(max(chg(idx),-limit),limit);
if (verbose && any(ochg ~= chg)),
disp(["Change in parameter(s): ", ...
sprintf("%d ", find(ochg ~= chg)), "were constrained"]);
endif
aprec=abs(pprec.*pbest); %---
%% ss=scalar sum of squares=sum((wt.*(y-f))^2).
if (any(abs(chg) > 0.1*aprec)),%--- % only worth evaluating function if
p=chg+pprev; % there is some non-miniscule change
%% apply bounds; preserving the direction of the attempted step
%% might lead to fixing _all_ parameters at their current value,
%% so decided not to do that and to simply reset parameters to
%% bounds
lidx = p < bounds(:, 1);
p(lidx) = bounds(lidx, 1);
uidx = p > bounds(:, 2);
p(uidx) = bounds(uidx, 2);
if (verbose && any (idx = lidx | uidx))
printf ("bounds apply for parameters number %s %i\n", \
sprintf ("%i, ", find (idx)(1:end - 1)), \
find (idx)(end));
endif
%%
f=feval(F,x,p);
r=wt.*(y-f);
ss=r'*r;
if ss<sbest,
pbest=p;
fbest=f;
sbest=ss;
endif
if ss<=sgoal,
break;
endif
endif %---
endfor
epsLlast = epsL;
if (verbose),
eval(plotcmd);
endif
if ss<eps,
break;
endif
aprec=abs(pprec.*pbest);
%% [aprec, chg, chgprev]
if (all(abs(chg) < aprec) && all(abs(chgprev) < aprec)),
cvg=1;
if (verbose),
fprintf("Parameter changes converged to specified precision\n");
endif
break;
else
chgprev=chg;
endif
if ss>sgoal,
break;
endif
endfor

%% set return values
p=pbest;
f=fbest;
ss=sbest;
cvg=((sbest>sgoal) || (sbest<=eps) || cvg);
if cvg ~= 1 , disp(" CONVERGENCE NOT ACHIEVED! "), endif

%% CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS
%% re-evaluate the Jacobian at optimal values
jac=feval(dFdp,x,f,p,dp,F,bounds);
msk = dp ~= 0;
n = sum(msk); % reduce n to equal number of estimated parameters
jac = jac(:, msk); % use only fitted parameters

%% following section is Ray Muzic's estimate for covariance and correlation
%% assuming covariance of data is a diagonal matrix proportional to
%% diag(1/wt.^2).
%% cov matrix of data est. from Bard Eq. 7-5-13, and Row 1 Table 5.1

if exist('sparse') % save memory
Q=sparse(1:m,1:m,1./wt.^2);
Qinv=sparse(1:m,1:m,wt.^2);
else
Q=diag((0*wt+1)./(wt.^2));
Qinv=diag(wt.*wt);
endif
resid=y-f; %un-weighted residuals
covr=resid'*Qinv*resid*Q/(m-n); %covariance of residuals
Vy=1/(1-n/m)*covr; % Eq. 7-13-22, Bard %covariance of the data

jtgjinv=inv(jac'*Qinv*jac); %argument of inv may be singular
covp=jtgjinv*jac'*Qinv*Vy*Qinv*jac*jtgjinv; % Eq. 7-5-13, Bard %cov of parm est
d=sqrt(abs(diag(covp)));
corp=covp./(d*d');

if exist('sparse')
covr=spdiags(covr,0);
stdresid=resid./sqrt(spdiags(Vy,0));
else
covr=diag(covr); % convert returned values to compact storage
stdresid=resid./sqrt(diag(Vy)); % compute then convert for compact storage
endif
Z=((m-n)*jac'*Qinv*jac)/(n*resid'*Qinv*resid);

%%% alt. est. of cov. mat. of parm.:(Delforge, Circulation, 82:1494-1504, 1990
%%disp('Alternate estimate of cov. of param. est.')
%%acovp=resid'*Qinv*resid/(m-n)*jtgjinv

%%Calculate R^2 (Ref Draper & Smith p.46)
r=corrcoef([y(:),f(:)]);
r2=r(1,2).^2;

%% if someone has asked for it, let them have it
if (verbose),
eval(plotcmd);
disp(" Least Squares Estimates of Parameters")
disp(p')
disp(" Correlation matrix of parameters estimated")
disp(corp)
disp(" Covariance matrix of Residuals")
disp(covr)
disp(" Correlation Coefficient R^2")
disp(r2)
sprintf(" 95%% conf region: F(0.05)(%.0f,%.0f)>= delta_pvec'*Z*delta_pvec",n,m-n)
Z
%% runs test according to Bard. p 201.
n1 = sum((f-y) < 0);
n2 = sum((f-y) > 0);
nrun=sum(abs(diff((f-y)<0)))+1;
if ((n1>10)&(n2>10)), % sufficent data for test?
zed=(nrun-(2*n1*n2/(n1+n2)+1)+0.5)/(2*n1*n2*(2*n1*n2-n1-n2)...
/((n1+n2)^2*(n1+n2-1)));
if (zed < 0),
prob = erfc(-zed/sqrt(2))/2*100;
disp([num2str(prob),"% chance of fewer than ",num2str(nrun)," runs."]);
else
prob = erfc(zed/sqrt(2))/2*100;
disp([num2str(prob),"% chance of greater than ",num2str(nrun)," runs."]);
endif
endif
endif

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by roessli bertrand :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

Hi,

Thanks, I will make the changes. It will take some time tough...

Regarding the "bounds" in nmsmax I did it in a primitive way
that is enough for my needs:

[f(1)] = dirn*feval(fun,x,varargin{:});

--> [f(1), param] = dirn*feval(fun,x,varargin{:});

where the function fun renormalise the parameters within 'bounds'
and send them back to nmsmax.

In my case these are not 'bounds' (b_min<=p<=b_max) ) in the
strict sense but I need that \sum{some_parameters} = 1.
(I can probably try something like 1-\sum_2{param}=param(1) with 0<=param(2,:)<=1)
I have to think a bit...

Regards,

Bertrand

---
http://broessli.redbubble.com/

--- On Tue, 8/4/09, Olaf Till <olaf.till@...> wrote:

From: Olaf Till <olaf.till@...>
Subject: Re: [OctDev] package contribution
To: octave-dev@...
Date: Tuesday, August 4, 2009, 1:02 PM

On Mon, Aug 03, 2009 at 01:49:45PM +0200, Olaf Till wrote:

> On Mon, Aug 03, 2009 at 01:07:23PM +0200, Olaf Till wrote:
> > I don't know if there are theoretical objections against bounds in
> > Levenberg-Marquardt optimization, but I also ended up introducing
> > bounds into leasqr for my use; I would think this is general
> > purpose.
>
> Sorry, I remember now that there were some problems. Often, a bound is
> desirable because some parameters do not produce useful results (of
> the user function) outside a certain range (e.g. a division by zero
> might occur). In such a case, the bounds should already be respected
> during gradient determination, not only later in the 'step'. This
> would mean that if a parameter is already at its bound, only a
> _one_-sided gradient can be computed for this parameter, but only on
> the 'side' away from the bound. Also, in two-sided bounds, the
> possible parameter-change for gradient estimation could become small
> (too small for some problems?). Since currently gradient estimation in
> leasqr is done by different function (dfdp), which knows nothing on
> bounds, implementation is not so easy. One could pass the
> bounds-argument also to dfdp, but there might exist user-replacements
> for dftp (there is an argument to leasqr just to choose this function)
> which do not honour the 'bounds' argument ... so people can not expect
> bounds to work without changing their dfdp-replacement.

Some new suggestions.

-- "bounds" in leasqr

Attached are new versions of leasqr.m and dfdp.m. They handle "bounds"
in steps and in gradient estimation. The "options" argument is used
for this and is now a structure of options. For backwards
compatibility, the old options-matrix is still recognized. A lot of
cosmetic changes have been made to make leasqr.m better editable under
emacs, so almost each line is changed and its not much use making
diffs :-(. I have tested it a bit. Since leasqr does not seem to have
an "assigned" active maintainer, I would commit it and change the
version number of the optim package, but Soeren please decide on that.

Roessli, if there is nothing wrong with the new leasqr version I
suggest you rewrite your code to make use of it, and make it depend on
the new version of the optim package. The computations according to
Bragg you had put inside it can surely be made _after_ calling leasqr?

-- "bounds" in nmsmax

I am not familiar with the method in nmsmax, but at first glance I do
not see where you introduced "bounds" ...

-- keeping parameters "fixed"

The way you have done this in nmsmax.m could break existing code since
the new argument is not the last. And it can not be made the last due
to the already existing varargin argument. (BTW, why not using an
argument which is a logical vector, true for each parameter to be
fixed?) To be on the safe side, I would suggest you to solve the
"fixing" problem with a wrapper, which could be left in your package
and need not necessarily put into the optim package. Something like
the following sketch (may have overlooked some details):

function [outp1, outp2, ...] = optim_wrapper (fixed, varargin)
## according to index "fixed", remove fixed parameters and possibly
## other related vector-elements
...
## replace name of (or handle to) user-function with an anonymous
## function like that (calls the here defined subfunction)
@ (x, p) local_func (x, p, user_func, fixed, values_of_fixed_parameters)
...
## call optimization routine
[oup1, outp2, ...] = call_optimizer (varargin(:));
## reintroduce fixed parameters into output arguments if necessary
...
## define subfunction
function ret = local_func (x, p, user_func, fixed, values_fixed)
## reintroduce fixed parameters into "p" according to "fixed" and
## "values_fixed"
...
## call user function
ret = feval (user_func, x, p);

-----Inline Attachment Follows-----

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july

-----Inline Attachment Follows-----

_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by Olaf Till :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

On Tue, Aug 04, 2009 at 04:47:15AM -0700, roessli bertrand wrote:
> Hi,
>
> Thanks, I will make the changes. It will take some time tough...

But please wait until the new leasqr.m is accepted before you change
the call to leasqr; I can not simply commit it but must wait for
approval :-)

Olaf

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by Søren Hauberg :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

tir, 04 08 2009 kl. 13:02 +0200, skrev Olaf Till:
> Attached are new versions of leasqr.m and dfdp.m. They handle "bounds"
> in steps and in gradient estimation.

Could you send diffs instead of entire files? That way it is much more
easy to see what kind of changes you have made.

> The "options" argument is used
> for this and is now a structure of options. For backwards
> compatibility, the old options-matrix is still recognized.

Sounds reasonable.

> A lot of
> cosmetic changes have been made to make leasqr.m better editable under
> emacs, so almost each line is changed and its not much use making
> diffs :-(.

Then could you please start out by doing a diff where you do not change
coding style. Once we have the content in place, you can always change
the coding style.

> I have tested it a bit. Since leasqr does not seem to have
> an "assigned" active maintainer, I would commit it and change the
> version number of the optim package, but Soeren please decide on that.

I know my basic optimisation, but I doubt that I'm the best candidate
for this job. So, as always, I appreciate comments from other people
(that way, we might actually end up making informed decisions).

Søren

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by Olaf Till :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

On Tue, Aug 04, 2009 at 09:57:28PM +0200, Søren Hauberg wrote:
> tir, 04 08 2009 kl. 13:02 +0200, skrev Olaf Till:
> > A lot of
> > cosmetic changes have been made to make leasqr.m better editable under
> > emacs, so almost each line is changed and its not much use making
> > diffs :-(.
>
> Then could you please start out by doing a diff where you do not change
> coding style. Once we have the content in place, you can always change
> the coding style.

Ok, attached is an SVN diff inside main/optim/inst/.

Olaf

Index: leasqr.m =================================================================== --- leasqr.m (revision 6076) +++ leasqr.m (working copy) @@ -46,17 +46,28 @@ % stol = scalar tolerance on fractional improvement in scalar sum of % squares = sum((wt.*(y-f))^2); default stol = .0001; % niter = scalar maximum number of iterations; default = 20; -% options = matrix of n rows (same number of rows as pin) containing -% column 1: desired fractional precision in parameter estimates. -% Iterations are terminated if change in parameter vector (chg) on two -% consecutive iterations is less than their corresponding elements -% in options(:,1). [ie. all(abs(chg*current parm est) < options(:,1)) -% on two consecutive iterations.], default = zeros(). -% column 2: maximum fractional step change in parameter vector. -% Fractional change in elements of parameter vector is constrained to be -% at most options(:,2) between sucessive iterations. -% [ie. abs(chg(i))=abs(min([chg(i) options(i,2)*current param estimate])).], -% default = Inf*ones(). +% options = structure, currently recognized fields are 'fract_prec', +% 'max_fract_change', and 'bounds'. For backwards compatibility, +% 'options' can also be a matrix whose first and second column +% contains the values of 'fract_prec' and 'max_fract_change', +% respectively. +% Field 'options.fract_prec': column vector (same length as 'pin') of +% desired fractional precisions in parameter estimates. Iterations +% are terminated if change in parameter vector (chg) on two +% consecutive iterations is less than their corresponding elements in +% 'options.fract_prec'. [ie. all(abs(chg*current parm est) < +% options.fract_prec) on two consecutive iterations.], default = +% zeros(). +% Field 'options.max_fract_change': column vector (same length as +% 'pin) of maximum fractional step changes in parameter vector. +% Fractional change in elements of parameter vector is constrained to +% be at most 'options.max_fract_change' between sucessive iterations. +% [ie. abs(chg(i))=abs(min([chg(i) +% options.max_fract_change(i)*current param estimate])).], default = +% Inf*ones(). +% Field 'options.bounds': two-column-matrix, one row for each +% parameter in 'pin'. Each row contains a minimal and maximal value +% for each parameter. Default: [-Inf, Inf] in each row. % % OUTPUT VARIABLES % f = column vector of values computed: f = F(x,p). @@ -122,6 +133,7 @@ % MATLAB version >= 4.0) in computation of covp, corp, stdresid. % Modified by Francesco Potortì % for use in Octave +% Added bounds feature to this file and to dfdp.m, Olaf Till 4-Aug-2009 % % References: % Bard, Nonlinear Parameter Estimation, Academic Press, 1974. @@ -158,27 +170,64 @@ y=y(:); wt=wt(:); pin=pin(:); dp=dp(:); %change all vectors to columns % check data vectors- same length? -m=length(y); n=length(pin); p=pin;[m1,m2]=size(x); +m=length(y); n=length(pin); [m1,m2]=size(x); if m1~=m ,error('input(x)/output(y) data must have same number of rows ') ,end; -if (nargin <= 9), - options=[zeros(n,1), Inf*ones(n,1)]; - nor = n; noc = 2; -else - [nor, noc]=size(options); - if (nor ~= n), - error('options and parameter matrices must have same number of rows'), - end; - if (noc ~= 2), - options=[options(:,1), Inf*ones(nor,1)]; - end; -end; -pprec=options(:,1); -maxstep=options(:,2); -% + %% processing of 'options' + pprec = zeros (n, 1); + maxstep = Inf * ones (n, 1); + bounds = cat (2, -Inf * ones (n, 1), Inf * ones (n, 1)); + if (nargin > 9) + if (ismatrix (options)) % backwards compatibility + tp = options; + options = struct ("fract_prec", tp(:, 1)); + if (columns (tp) > 1) + options.max_fract_change = tp(:, 2); + endif + endif + if (isfield (options, "fract_prec")) + pprec = options.fract_prec; + if (rows (pprec) != n || columns (pprec) != 1) + error ("fractional precisions: wrong dimensions"); + endif + endif + if (isfield (options, "max_fract_change")) + maxstep = options.max_fract_change; + if (rows (maxstep) != n || columns (maxstep) != 1) + error ("maximum fractional step changes: wrong dimensions"); + endif + endif + if (isfield (options, "bounds")) + bounds = options.bounds; + if (rows (bounds) != n || columns (bounds) != 2) + error ("bounds: wrong dimensions"); + endif + idx = bounds(:, 1) > bounds(:, 2); + tp = bounds(idx, 2); + bounds(idx, 2) = bounds(idx, 1); + bounds(idx, 1) = tp; + if (any (idx = bounds(:, 1) == bounds(:, 2))) + warning ("lower and upper bounds identical for some parameters, setting the respective elements of 'dp' to zero"); + dp(idx) = 0; + endif + if (any (idx = pin < bounds(:, 1))) + warning ("some initial parameters set to lower bound"); + pin(idx) = bounds(idx, 1); + endif + if (any (idx = pin > bounds(:, 2))) + warning ("some initial parameters set to upper bound"); + pin(idx) = bounds(idx, 2); + endif + endif + endif + if (all (dp == 0)) + error ("no free parameters"); + endif + % set up for iterations % +p = pin; f=feval(F,x,p); fbest=f; pbest=p; r=wt.*(y-f); ss=r'*r; @@ -193,7 +242,7 @@ % for iter=1:niter, pprev=pbest; - prt=feval(dFdp,x,fbest,pprev,dp,F); + prt=feval(dFdp,x,fbest,pprev,dp,F,bounds); r=wt.*(y-fbest); sprev=sbest; sgoal=(1-stol)*sprev; @@ -234,6 +283,20 @@ % ss=scalar sum of squares=sum((wt.*(y-f))^2). if (any(abs(chg) > 0.1*aprec)),%--- % only worth evaluating function if p=chg+pprev; % there is some non-miniscule change + %% apply bounds; preserving the direction of the attempted step + %% might lead to fixing _all_ parameters at their current value, + %% so decided not to do that and to simply reset parameters to + %% bounds + lidx = p < bounds(:, 1); + p(lidx) = bounds(lidx, 1); + uidx = p > bounds(:, 2); + p(uidx) = bounds(uidx, 2); + if (verbose && any (idx = lidx | uidx)) + printf ("bounds apply for parameters number %s %i\n", \ + sprintf ("%i, ", find (idx)(1:end - 1)), \ + find (idx)(end)); + endif + %% f=feval(F,x,p); r=wt.*(y-f); ss=r'*r; @@ -280,7 +343,7 @@ % CALC VARIANCE COV MATRIX AND CORRELATION MATRIX OF PARAMETERS % re-evaluate the Jacobian at optimal values -jac=feval(dFdp,x,f,p,dp,F); +jac=feval(dFdp,x,f,p,dp,F,bounds); msk = dp ~= 0; n = sum(msk); % reduce n to equal number of estimated parameters jac = jac(:, msk); % use only fitted parameters Index: dfdp.m =================================================================== --- dfdp.m (revision 6076) +++ dfdp.m (working copy) @@ -15,7 +15,7 @@ % You should have received a copy of the GNU General Public License % along with this program; If not, see <http://www.gnu.org/licenses/>. -function prt=dfdp(x,f,p,dp,func) +function prt=dfdp(x,f,p,dp,func,bounds) % numerical partial derivatives (Jacobian) df/dp for use with leasqr % --------INPUT VARIABLES--------- % x=vec or matrix of indep var(used as arg to func) x=[x0 x1 ....] @@ -27,26 +27,55 @@ % dp(j)=0 sets corresponding partials to zero; i.e. holds p(j) fixed % func=string naming the function (.m) file % e.g. to calc Jacobian for function expsum prt=dfdp(x,f,p,dp,'expsum') +% bounds=two-column-matrix of lower and upper bounds for parameters %----------OUTPUT VARIABLES------- % prt= Jacobian Matrix prt(i,j)=df(i)/dp(j) %================================ m=size(x,1); if (m==1), m=size(x,2); end %# PAK: in case #cols > #rows n=length(p); %dimensions -ps=p; prt=zeros(m,n);del=zeros(n,1); % initialise Jacobian to Zero -for j=1:n - del(j)=dp(j) .*p(j); %cal delx=fract(dp)*param value(p) - if p(j)==0 - del(j)=dp(j); %if param=0 delx=fraction - end - p(j)=ps(j) + del(j); - if del(j)~=0, f1=feval(func,x,p); - if dp(j) < 0, prt(:,j)=(f1-f)./del(j); - else - p(j)=ps(j)- del(j); - prt(:,j)=(f1-feval(func,x,p))./(2 .*del(j)); - end - end - p(j)=ps(j); %restore p(j) -end -return +prt=zeros(m,n);del=zeros(n,1); % initialise Jacobian to Zero +del = dp .* p; %cal delx=fract(dp)*param value(p) +idx = p == 0; +del(idx) = dp(idx); %if param=0 delx=fraction +idx = dp > 0; +del(idx) = abs (del(idx)); % not for one-sided intervals, changed + % direction of intervals could change + % behavior of optimization without bounds +min_del = min (abs (del), bounds(:, 2) - bounds(:, 1)); + for j=1:n + ps = p; + if (dp(j)~=0) + if (dp(j) < 0) + ps(j) = p(j) + del(j); + if (ps(j) < bounds(j, 1) || ps(j) > bounds(j, 2)) + t_del1 = max (bounds(j, 1) - p(j), - abs (del(j))); % + %non-positive + t_del2 = min (bounds(j, 2) - p(j), abs (del(j))); % + %non-negative + if (- t_del1 > t_del2) + del(j) = t_del1; + else + del(j) = t_del2; + endif + ps(j) = p(j) + del(j); + endif + prt(:, j) = (feval (func, x, ps) - f) / del(j); + else + if (p(j) - del(j) < bounds(j, 1)) + tp = bounds(j, 1); + ps(j) = tp + min_del(j); + elseif (p(j) + del(j) > bounds(j, 2)) + ps(j) = bounds(j, 2); + tp = ps(j) - min_del(j); + else + ps(j) = p(j) + del(j); + tp = p(j) - del(j); + min_del(j) = 2 * del(j); + endif + f1 = feval (func, x, ps); + ps(j) = tp; + prt(:, j) = (f1 - feval (func, x, ps)) / (min_del(j)); + endif + endif + endfor

Re: package contribution

by roessli bertrand :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

Hello,

Just have a comment which is also a question concerning dfdp.m

The help states
"dp= fractional increment of p for numerical derivatives"
which is fine, but the problem comes that usually nobody knows
which step to use to calculate the derivatives.

Would it no be better to choose the step in a similar way than it is done
in __bgsmin.cc?:

SQRT_EPS = sqrt(DBL_EPSILON);
        diff = exp(log(DBL_EPSILON)/3.0);
         // get 1st derivative by central difference
        for (j=0; j<k; j++) {
                p = parameter(j);
                // determine delta for finite differencing
                test = (fabs(p) + SQRT_EPS) * SQRT_EPS > diff;
                if (test) delta = (fabs(p) + SQRT_EPS) * SQRT_EPS;
                else delta = diff;
                // right side
                parameter(j) = p + delta;
                success = __bfgsmin_obj(obj_right, f, f_args, parameter, minarg);
                if (!success) error("__numgradient: objective function failed, can't compute numeric gradient");
                // left size
                parameter(j) = p - delta;
                success = __bfgsmin_obj(obj_left, f, f_args, parameter, minarg);
                if (!success) error("__numgradient: objective function failed, can't compute numeric gradient");                parameter(j) = p; // restore original parameter for next round
                g(j) = (obj_right - obj_left) / (2.0*delta);
        }

Best regards,

Bertrand
---
http://broessli.redbubble.com/

--- On Tue, 8/4/09, Olaf Till <olaf.till@...> wrote:

From: Olaf Till <olaf.till@...>
Subject: Re: [OctDev] package contribution
To: octave-dev@...
Date: Tuesday, August 4, 2009, 1:02 PM

On Mon, Aug 03, 2009 at 01:49:45PM +0200, Olaf Till wrote:

> On Mon, Aug 03, 2009 at 01:07:23PM +0200, Olaf Till wrote:
> > I don't know if there are theoretical objections against bounds in
> > Levenberg-Marquardt optimization, but I also ended up introducing
> > bounds into leasqr for my use; I would think this is general
> > purpose.
>
> Sorry, I remember now that there were some problems. Often, a bound is
> desirable because some parameters do not produce useful results (of
> the user function) outside a certain range (e.g. a division by zero
> might occur). In such a case, the bounds should already be respected
> during gradient determination, not only later in the 'step'. This
> would mean that if a parameter is already at its bound, only a
> _one_-sided gradient can be computed for this parameter, but only on
> the 'side' away from the bound. Also, in two-sided bounds, the
> possible parameter-change for gradient estimation could become small
> (too small for some problems?). Since currently gradient estimation in
> leasqr is done by different function (dfdp), which knows nothing on
> bounds, implementation is not so easy. One could pass the
> bounds-argument also to dfdp, but there might exist user-replacements
> for dftp (there is an argument to leasqr just to choose this function)
> which do not honour the 'bounds' argument ... so people can not expect
> bounds to work without changing their dfdp-replacement.

Some new suggestions.

-- "bounds" in leasqr

Attached are new versions of leasqr.m and dfdp.m. They handle "bounds"
in steps and in gradient estimation. The "options" argument is used
for this and is now a structure of options. For backwards
compatibility, the old options-matrix is still recognized. A lot of
cosmetic changes have been made to make leasqr.m better editable under
emacs, so almost each line is changed and its not much use making
diffs :-(. I have tested it a bit. Since leasqr does not seem to have
an "assigned" active maintainer, I would commit it and change the
version number of the optim package, but Soeren please decide on that.

Roessli, if there is nothing wrong with the new leasqr version I
suggest you rewrite your code to make use of it, and make it depend on
the new version of the optim package. The computations according to
Bragg you had put inside it can surely be made _after_ calling leasqr?

-- "bounds" in nmsmax

I am not familiar with the method in nmsmax, but at first glance I do
not see where you introduced "bounds" ...

-- keeping parameters "fixed"

The way you have done this in nmsmax.m could break existing code since
the new argument is not the last. And it can not be made the last due
to the already existing varargin argument. (BTW, why not using an
argument which is a logical vector, true for each parameter to be
fixed?) To be on the safe side, I would suggest you to solve the
"fixing" problem with a wrapper, which could be left in your package
and need not necessarily put into the optim package. Something like
the following sketch (may have overlooked some details):

function [outp1, outp2, ...] = optim_wrapper (fixed, varargin)
## according to index "fixed", remove fixed parameters and possibly
## other related vector-elements
...
## replace name of (or handle to) user-function with an anonymous
## function like that (calls the here defined subfunction)
@ (x, p) local_func (x, p, user_func, fixed, values_of_fixed_parameters)
...
## call optimization routine
[oup1, outp2, ...] = call_optimizer (varargin(:));
## reintroduce fixed parameters into output arguments if necessary
...
## define subfunction
function ret = local_func (x, p, user_func, fixed, values_fixed)
## reintroduce fixed parameters into "p" according to "fixed" and
## "values_fixed"
...
## call user function
ret = feval (user_func, x, p);

-----Inline Attachment Follows-----

------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july

-----Inline Attachment Follows-----

_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Re: package contribution

by Olaf Till :: Rate this Message:

Reply to Author | View Threaded | Show Only this Message

On Thu, Aug 06, 2009 at 03:37:25AM -0700, roessli bertrand wrote:

> Hello,
>
> Just have a comment which is also a question concerning dfdp.m
>
> The help states
> "dp= fractional increment of p for numerical derivatives"
> which is fine, but the problem comes that usually nobody knows
> which step to use to calculate the derivatives.
>
> Would it no be better to choose the step in a similar way than it is done
> in __bgsmin.cc?:
>
> SQRT_EPS = sqrt(DBL_EPSILON);
>         diff = exp(log(DBL_EPSILON)/3.0);
>          // get 1st derivative by central difference
>         for (j=0; j<k; j++) {
>                 p = parameter(j);
>                 // determine delta for finite differencing
>                 test = (fabs(p) + SQRT_EPS) * SQRT_EPS > diff;
>                 if (test) delta = (fabs(p) + SQRT_EPS) * SQRT_EPS;
>                 else delta = diff;
>                 // right side
>                 parameter(j) = p + delta;
>                 success = __bfgsmin_obj(obj_right, f, f_args, parameter, minarg);
>                 if (!success) error("__numgradient: objective function failed, can't compute numeric gradient");
>                 // left size
>                 parameter(j) = p - delta;
>                 success = __bfgsmin_obj(obj_left, f, f_args, parameter, minarg);
>                 if (!success) error("__numgradient: objective function failed, can't compute numeric gradient");                parameter(j) = p; // restore original parameter for next round
>                 g(j) = (obj_right - obj_left) / (2.0*delta);
>         }

dfdp.m is meant to be called by leasqr.m, which chooses a default of
.001 for dp. This may not be optimal in some cases, but why do you
think the result of the code above is better?

Olaf
------------------------------------------------------------------------------
Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day
trial. Simplify your report design, integration and deployment - and focus on
what you do best, core application coding. Discover what's new with
Crystal Reports now. http://p.sf.net/sfu/bobj-july
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

< Prev | 1 - 2 | Next >