Contents

Diagnostic Results

This is a script to show the results of a run for IGUG It is called directly from IGUG but if you wish to run with given results you can upload a .mat file of data To upload data comment out and give the filename for results generated from IGUG load(Resultsfilename); filename_results=Resultsfilename; September 16, 2018, Rosemary Renaut December 19, 2018. Modified to adjust the range for data fitting. Also the mass estimate is plotted over time.

load('DataN410K1000_1')

Should the results be saved as figures and spreadsheet with data fitting

button = questdlg('Do you want to do linear regression analysis for diagnostics?','Diagnostics','Yes','No','Yes');
datafit=0;savefigs=0;spreadsheet=0;datafitrange=0;
if (exist('EData','var') | (~isempty('EData'))) realdata=0; else realdata=1;end
switch button, case "Yes", datafit=1;

    button = questdlg('Do you want to give the range for the data fitting?','Data Fit Range','Yes','No','Yes');

    switch button, case "Yes", datafitrange=1;end
end
button = questdlg('Do you want to save figures','Diagnostics','Yes','No','Yes');
switch button, case "Yes", savefigs=1;end

button = questdlg('Do you want to save results to spreadsheet','Diagnostics','Yes','No','Yes');
switch button, case "Yes", spreadsheet=1;end

Plot the results and do datafitting if requested

BestTetaEnd=[];BestMisfitEnd=[];BestFitEnd=[];startk=1;
ll=length(lambdai);
if (ll <= 4 ), px=1;py=ll; elseif ( ll==5 |  ll==6), px=3;py=2; elseif (ll <=9), px=3;py=3; else px=4;py=3;end
for ii=1:length(lambdai)
    % Info for the table
    BestTetaEnd=[BestTetaEnd;BestTetaAll(ii,count(ii))];
    BestMisfitEnd=[BestMisfitEnd;BestMisfitAll(ii,count(ii))];
    BestFitEnd=[BestFitEnd;BestFitAll(ii,count(ii))];
    k=count(ii)+1;
    % End of info for the table
    if datafit
        krange=[(min(75,floor((count(ii)+1)/2))),(count(ii))];
        if datafitrange
            button = questdlg(['Do you want to pick the range for the data fit for lambda = ',num2str(lambdai(ii))],'Data fit for lambda','Yes','No','Yes');

            switch button
                case "Yes"
                    kr={'75','200'};
                    prompt={'Lower value for k','Upper value for k'};
                    kr=inputdlg(prompt,'[kstart, kend]',[1,50],kr);
                    ks=str2num(cell2mat(kr(1)));ke=str2num(cell2mat(kr(2)));
                    krange=[ks,ke];krange=sort(krange,'ascend');

            end
        end
        startk=krange(1);endk=krange(2);if endk>k, endk=k-1;end
        [datafit1(:,ii),datafit2(:,ii),datafit3(:,ii),res1(ii),res2(ii),res3(ii),yhat1,yhat2,yhat3]=...
            Regression_Analysis(BestFitAll(ii,1:k-1), (BestTetaAll(ii,1:k-1)),BestMisfitAll(ii,1:k-1),endk,startk);
        figure(4),subplot(px,py,ii), plot(1:k-1,log(BestTetaAll(ii,1:k-1)),'--kd',startk:endk-1,yhat1,'*b','LineWidth',0.5);
        legend('Results','DataFit');title(['\lambda = ',num2str(lambdai(ii)),' R2 = ',num2str(res1(ii))]);
        figure(5),subplot(px,py,ii), plot(1:k-1,BestMisfitAll(ii,1:k-1),'--kd',startk:endk-1,yhat2,'*b','LineWidth',0.5);
        legend('Results','DataFit');title(['\lambda = ',num2str(lambdai(ii)),' R2 = ',num2str(res2(ii))])
        figure(6),subplot(px,py,ii), plot(1:k-1,BestFitAll(ii,1:k-1),'--kd',startk:endk-1,yhat3,'*b','LineWidth',0.5);
        legend('Results','DataFit');title(['\lambda = ',num2str(lambdai(ii)),' R2 = ',num2str(res3(ii))])
        figure(7),subplot(px,py,ii),contourf(Estat,Nstat,(reshape(gz(ii,:),nse,nsn))');shading flat;
        title(['\lambda = ',num2str(lambdai(ii))]);
    else
        figure(4),subplot(px,py,ii), plot(1:k-1,log(BestTetaAll(ii,1:k-1)),'LineWidth',0.5);
        legend('Results' );title(['\lambda = ',num2str(lambdai(ii))]);
        figure(5),subplot(px,py,ii), plot(1:k-1,BestMisfitAll(ii,1:k-1),'LineWidth',0.5);
        legend('Results' );title(['\lambda = ',num2str(lambdai(ii))])
        figure(6),subplot(px,py,ii), plot(1:k-1,BestFitAll(ii,1:k-1),'LineWidth',0.5);
        legend('Results' );title(['\lambda = ',num2str(lambdai(ii))])
        figure(7),subplot(px,py,ii),contourf(Estat,Nstat,(reshape(gz(ii,:),nse,nsn))');shading flat;
        title(['\lambda = ',num2str(lambdai(ii))]);
    end
end
figure(4),suptitle(['\fontsize{16} \bf \Theta(p)']);set(gcf,'Units','Centimeter','Position',[1 5 25 25]);
figure(5),suptitle(['\fontsize{16} \bf \Phi(p)']);set(gcf,'Units','Centimeter','Position',[26 5 25 25]);
figure(6),suptitle(['\fontsize{16} \bf \Gamma(p)']);set(gcf,'Units','Centimeter','Position',[51 5 25 25]);
figure(7),suptitle(['\fontsize{16} \bf Predicted Data']);set(gcf,'Units','Centimeter','Position',[76 5 25 25]);

Plot the Mass Point Distribution

figure(3),clf
for ii=1:length(lambdai)
    subplot(px,py,ii),
    PlotDomainStructure(Fac,Ver,East1,East2,North1,North2,Depth1,Depth2,prismdef,noc)
    ScatterPlotResults(FpAll(ii,:),M);
    if datafit
        title(['\lambda = ',num2str(lambdai(ii)),' R2 = ',num2str(res1(ii))])
    else
        title(['\lambda = ',num2str(lambdai(ii))])
    end
end
set(gcf,'Units','Centimeter','Position', [3 1 30 30]);

Plot the Mass with iteration

figure(9),clf
for ii=1:length(lambdai)
    subplot(px,py,ii),plot(1:count(ii),(BestMassAll(ii,1:count(ii))),'LineWidth',1.0);
     if datafit
        title(['\lambda = ',num2str(lambdai(ii)),' R2 = ',num2str(res1(ii))])
    else
        title(['\lambda = ',num2str(lambdai(ii))])
    end
end
suptitle(['\fontsize{16} \bf mass']);
set(gcf,'Units','Centimeter','Position', [33,45,35,10]);figprops

If figures should be saved then do printfig on each

if ~exist('filename_results')| (isempty('filename_results'))
    if realdata
        filename_results = ['ResultsRealM',int2str(M)];
    else
        filename_results = ['ResultsDataM',int2str(M)];
    end
end
if savefigs
    for i=[3:6,9]
        figure(i)
        printfig(['Figure',int2str(i),filename_results]);
    end
    if (exist('EData') && (~isempty('EData')))
        figure(7), printfig(['Figure',int2str(7),filename_results]);end
end

Show the results as a table dependent on whether datafit has been done

mm=N+sqrt(2*N);
if datafit
    T=table(M*ones(length(lambdai),1),lambdai(:),count(:),mass(:),BestTetaEnd(:),BestMisfitEnd(:),BestMisfitEnd(:)/mm,res1(:),time_lambda(:),...; %time_lambda(:),res1(:),...
        'VariableNames',{'M','lambda','kend','mass','Theta','Phi','Phi_over_chi','R2Theta','Time'})
%     T=table(mass(:),count(:),BestTetaEnd,BestMisfitEnd/mm,BestFitEnd,lambdai(:),datafit1',res1(:),datafit2',res2(:),datafit3',res3(:),...
%         'VariableNames',{'mass','kend','Theta','Phi_over_chi','Gamma','lambda',...
%         'logdatafitTheta','r2theta','lindatafitPhi','r2Phi','lindatafitGamma','r2Gamma'})
else
    T=table(mass(:),count(:),BestTetaEnd,BestMisfitEnd/mm,BestFitEnd,lambdai(:),...
        'VariableNames',{'mass','kend','Theta','Phi_over_chi','Gamma','lambda'})
end
T =

  12×9 table

    M     lambda    kend       mass         Theta        Phi      Phi_over_chi    R2Theta      Time 
    __    ______    ____    __________    __________    ______    ____________    ________    ______

    10       200    1000    1.0349e+11       0.47455     88380       66.883        0.19935    429.76
    10        20    1000    1.1264e+11        32.827     48462       36.674        0.29108    430.96
    10         2    1000    1.1516e+11        5.4948    2780.7       2.1043        0.12168    436.92
    10         1    1000    1.2035e+11         20.66    3073.1       2.3256        0.14792    449.65
    10       0.5    1000    1.1829e+11          12.6    2210.1       1.6725        0.47996    462.33
    10       0.2    1000    1.1326e+11        9.4254    1735.4       1.3133        0.88649     452.4
    10       0.1    1000    1.1159e+11        29.882    1581.2       1.1966        0.48692    444.59
    10      0.05    1000    1.0978e+11         57.29    1370.4       1.0371        0.55778    431.45
    10      0.02    1000    1.1005e+11         471.1    1392.8        1.054       0.072785     429.6
    10     0.002    1000    1.1142e+11        3469.7    1368.8       1.0359        0.46337     436.5
    10    0.0002    1000    1.1174e+11         20837    1380.2       1.0445        0.51551    443.81
    10     2e-05    1000    1.1667e+11    1.3195e+05    1444.8       1.0934        0.33193    450.38

if spreadsheet
    writetable(T,['Table',filename_results,'.xls']);
end