%First pass at reviving the code to optimize mask-analysis %code format compact clear all %test data from 1978 paper: %scattering matrix smat = ... [ 392 33 81 51 3 2 0 1 1 39 91 112 59 12 8 1 0 1 61 64 71 157 25 14 1 2 3 4 10 11 19 39 22 6 5 4 3.6 .7 1.4 4.6 .7 3.3 3.1 5.9 3.5] ; smatRowNames = ... [ ' nuc', ' mit', ' mitrym', ... ' ocyt', ' sch', ' schrim', ' schax', ' axon', ' axonrim']; smatColNames = [' nuc', ' mit', ' ocyt', ' sch', ' axon']; % real grains observed = ... [ 11 21 19 24 12 17 7 9 9 ]; % print observed fprintf (1,'%s\n', 'Observed grain distribution:'); fprintf (1,'%s', smatRowNames); fprintf (1,'\n'); fprintf (1,'%8.2f', observed); fprintf (1,'\n\n') %print scatter matrix fprintf (1,'%s\n', 'Generated source-to-grain matrix (grain comp across)') for row=1:size(smat,1) fprintf (1,'%8.2f', smat(row,:)); fprintf (1,'\n'); end fprintf (1,'\n'); %print sums of cols and relative sizes fprintf (1,'%s\n', 'Sums of scattering matrix columns:'); fprintf (1,'%8.2f', sum(smat,1)); fprintf (1,'\n\n') fprintf (1,'%s\n', 'Relative grain comparment sizes (percent):'); fprintf (1,'%8.2f', 100*sum(smat,1)/sum(sum(smat)) ); fprintf (1,'\n\n') fprintf (1,'%s\n', 'Sum of all matrix elements:'); fprintf (1,'%8.2f', sum(sum(smat)) ); fprintf (1,'\n\n') %get initial estimate for source distribution and print it Dinitial = sum(observed) / sum(sum(smat)); fprintf (1,'%s\n', 'Average observed grain density:'); fprintf (1,'%10.4f', Dinitial ); fprintf (1,'\n\n'); %make an initial row-vector for the density, % of size equal to the number of rows in % smat, with every element equal to Dinitial. density00 = ( ones(1,size(smat,1))* Dinitial ); %make a initial row-vector for the density determined %from the actual observed grains and smat density0 = density00 % = observed * pinv(smat); fprintf (1,'%s\n', 'Initial density estimate:'); fprintf (1,'%s', smatColNames); fprintf (1,'\n'); fprintf (1,'%8.4f', density0 ); fprintf (1,'\n\n'); %function to be optimized is in chisq.m %First the unconstrained version options = []; fixdensity = ones(1,size(smat,1)); %don't constrain any variable ChiMin = fmins ('chisq', density0, options, [], smat, observed, fixdensity); fprintf (1,'%s\n', 'Optimized source density--no constraints:'); fprintf (1,'%s', smatColNames); fprintf (1,'\n'); fprintf (1,'%8.4f', ChiMin); fprintf (1,'\n'); model = ChiMin*smat; %get the source strength error estimate Nsource = size(smat,1); Mgrain = size(smat, 2); %get Q for i=1:Nsource for k=1:Nsource sumterm=0; for j=1:Mgrain sumterm = sumterm + smat(i,j)*smat(k,j)/model(j); end Q(i,k) = sumterm; end end Qinv = inv(Q); %get P for i=1:Nsource for j=1:Mgrain sumterm=0; for k=1:Nsource sumterm = sumterm + Qinv(i,k)*smat(k,j)/model(j); end P(i,j) = sumterm; end end %get the error for i=1:Nsource sumterm=0; for j=1:Mgrain sumterm = sumterm + P(i,j)^2 * model(j); end dmodel(i) = sqrt(sumterm); end fprintf (1,'%s\n', 'Computed grain distribution error ranges:'); fprintf (1,'%8.4f', dmodel); fprintf (1,'\n'); fprintf (1,'%s\n', 'Computed grain distribution from optimum sources:'); fprintf (1,'%s', smatRowNames); fprintf (1,'\n'); fprintf (1,'%8.2f', model); fprintf (1,'\n'); fprintf (1,'%s\n', 'Components of ChiSq from optimum sources:'); fprintf (1,'%8.2f',((observed - model) .^ 2) ./ model ); fprintf (1,'\n'); fprintf (1,'%s\n', 'Total ChiSq from optimum sources:'); fprintf (1,'%8.2f', sum ( ((observed - model) .^ 2) ./ model )); fprintf (1,'\n\n'); %Now the code to constrain to all densities > 0.0 count = 0; while ( sum(ChiMin < 0) > 0 & count<10 ) %then there must be at least one neg element fixdensity = (ChiMin > .001*density00); density0 = density0 .* fixdensity; ChiMin = fmins ('chisq', density0, options, [], smat, observed, fixdensity); count = count + 1; end ChiMin = ChiMin .* fixdensity; %Now force the zero components fprintf (1,'%s\n', 'Optimized source density--positive constraint:'); fprintf (1,'%s', smatColNames); fprintf (1,'\n'); fprintf (1,'%8.4f', ChiMin); fprintf (1,'\n') model = ChiMin*smat; %get the source strength error estimate Nsource = size(smat,1); Mgrain = size(smat, 2); %get Q for i=1:Nsource for k=1:Nsource sumterm=0; for j=1:Mgrain sumterm = sumterm + smat(i,j)*smat(k,j)/model(j); end Q(i,k) = sumterm; end end Qinv = inv(Q); %get P for i=1:Nsource for j=1:Mgrain sumterm=0; for k=1:Nsource sumterm = sumterm + Qinv(i,k)*smat(k,j)/model(j); end P(i,j) = sumterm; end end %get the error for i=1:Nsource sumterm=0; for j=1:Mgrain sumterm = sumterm + P(i,j)^2 * model(j); end dmodel(i) = sqrt(sumterm); end fprintf (1,'%s\n', 'Computed grain distribution error ranges:'); fprintf (1,'%8.4f', dmodel); fprintf (1,'\n'); fprintf (1,'%s\n', 'Computed grain distribution from optimum positive sources:'); fprintf (1,'%s', smatRowNames); fprintf (1,'\n'); fprintf (1,'%8.2f', model); fprintf (1,'\n'); fprintf (1,'%s\n', 'Components of ChiSq from optimum positive sources:'); fprintf (1,'%8.2f',((observed - model) .^ 2) ./ model ); fprintf (1,'\n'); fprintf (1,'%s\n', 'Total ChiSq from optimum positive sources:'); fprintf (1,'%8.2f', sum ( ((observed - model) .^ 2) ./ model )); fprintf (1,'\n\n');