monteTL_3D_5minbin_kogia.m

% Montecarlo  - 3D beam, kogia detectability model.
% Kait Frasier
% 7/23/2013
%{
% Modeling individual click detectability = "group counting method"
% This means that every point generated by the model just represents a
% group of animals. The parameters of each are drawn from a series of
% probability distributions, which ultimately combine to dictate a recieved
% level (RL), at which point a threshold is used to decide whether or not
% the group is "heard".

%}
clearvars
% if matlabpool('size') == 0
%     matlabpool
% else
%     matlabpool close force
%     matlabpool
% end
TLprofile_directory = 'E:\JAH\Kogia\kogia_TL_models\';
saveDir = 'E:\JAH\Kogia\Plots\';

spVec = {'kogia'};
siteVec = {'DT','GC','MC'};
for itSp = 1:length(spVec)
    % Check if output directory exits
    if ~exist(saveDir,'dir')
        mkdir(saveDir)
    end
    species = spVec{itSp};
    for itSite = 1:length(siteVec)
        site = siteVec{itSite};
        
        polarFile = fullfile(TLprofile_directory,sprintf('%skogia_Jan\\118kHz\\%szoom1_118_3DTL_2.mat',site,site));
        
        outdir = fileparts(polarFile);
        load(polarFile)
        
        if ~exist('botDepthSort','var')
            botDepthSort = botDepth_interp(IX,:);
        end
        for preIt = 1:length(sortedTLVec)
            sortedTLVec{preIt}(:,1) = sortedTLVec{preIt}(:,2);
        end
        % sort the bottom vectors by increasing angle, for some reason they
        % were not - file order issue, way back.
        
        n = 500; % the number of model runs that will feature in the probability distribution
        N = 10000; % simulate this many points per model run

        % variables to pick from a distribution for CV estimation
        diveDepth_mean = 700 + 100*rand(n,1); % mean dive altitude is somewhere between 175 and 225m
        diveDepth_std = 25 + 25*rand(n,1);% dive depth std. dev. is 10 to 20m
        SL_std = 2 + 3*rand(n,1); % add 1 to 3 db std dev to source level.
        
        
        % %
        SL_mean = 210 + 5*rand(n,1); % mean source level is between 205 and 210 dB pp - gervais
        clickStart_mean = 50 + 50*rand(n,1); % depth in meters at which clicking starts, from tyack 2006 for cuvier's;
        clickStart_std = 10 + 10*rand(n,1); % depth in meters at which clicking starts, from tyack 2006 for cuvier's;
        descAngle_std = 5 + 5*rand(n,1); % std deviation of descent angle in deg, from tyack 2006 for cuvier's;
        descAngle_mean = 72 + 5*rand(n,1);% mean descent angle in deg, from tyack 2006 for cuvier's;
        
        % directivity_mean = 30 + 5*rand(n,1); %beam directivity
        % directivity_std = 2 + 3*rand(n,1);% deviation of directivity is between 3 and 5dB
        minBeamAmp_mean = 26 + 3*rand(n,1); % minimum off-axis dBs down from peak
        minBeamAmp_std = 2+3*rand(n,1);
        
        numAngle = length(thisAngle);
        maxRange = 1000; % in meters
        thresh = 128;% click detection threshold (amplitude in dB pp)
        rr_int = round(rr(2)-rr(1)); % figure out what the range step size is
        nrr_new = rr_int*nrr;
        rr_new = 0:rr_int:nrr_new; % What are the real values of the range vector? (in m)
        pDetTotal = nan(n,1);
        
        binVec = 0:100:maxRange;
        binnedPercDet = nan(n,length(binVec)-1);
        binnedCounts = nan(n,length(binVec)-1);
        rotHorizDeg = 140 + 20*rand(n,1);
        rotHorizDeg_std = 10 + 10*rand(n,1);
        rotVertDegForage = 55 + 10*rand(n,1);
        rotVertDegForageStd = 5 + 5*rand(n,1);
        rotVertDegDive = 40 + 20*rand(n,1);
        rotVertDegDiveStd = 5 + 10*rand(n,1);
        
        descentPerc = .10 + .05*rand(n,1);
        
        % set up various depth limits to see how different distributions affect the
        % p(det) (ie. how sensitive is the model).
        % minDepth_vec = 200; %[0, 500, 800, 1000, 0]; in meters
        % maxDepth_vec = 1000; %[1000, 1500, 1200, 2000, 2000]; in meters
        
        
        %for itr0 = 1:length(minDepth_vec) % depth limit loop
        %    minDepth = minDepth_vec(itr0);
        %    maxDepth = maxDepth_vec(itr0);
        
        for itr_n = 1:n  % number of simulations loop
            if rem(itr_n,100) == 0
                fprintf('TL computation %d of %d\n', itr_n, n)
            end
            %%%%% Location Computation %%%%%
            % rand location
            randVec = ceil(rand(2,N)'.*repmat([2*maxRange, 2*maxRange], [N, 1]))...
                - repmat([maxRange, maxRange], [N, 1]);
            [theta, rho] = cart2pol(randVec(:,1),randVec(:,2));  % convert to polar coord.
            clear randVec  % trying to save on memory
            
            % trim out the locations that are beyond the max range (corners of the
            % 2*maxRange X 2*maxRange square, since now we are using a pi*maxRange^2
            % circle)
            jjj = 1;
            rho2 = [];
            theta2 = [];
            for iii = 1:length(rho)
                if rho(iii) < maxRange
                    rho2(jjj,1) = rho(iii);
                    theta2(jjj,1) = theta(iii);
                    jjj = jjj+1;
                end
            end
            thetaDeg = make360(theta2*180/pi);
            clear theta
            
            % go from angle to ref indices - pulled this into a function
            % because it happens a few times.
            [angleRef,radRef] = angle_ref_comp(thetaDeg,rho2,thisAngle);
            
            %%%%% Depth Computation %%%%%
            % Compute bottom depth at each randomly selected point
            count0 = 1;
            tempDepth = zeros(size(angleRef));
            keepPoint = ones(size(angleRef));
            
            diveDepthRef = diveDepth_mean(itr_n) + diveDepth_std(itr_n)...
                *randn(size(angleRef));
            % Find points that are above surface or below bottom and correct them
            flyingWhaleIdx = find(diveDepthRef>=sd | diveDepthRef<1);
            while ~isempty(flyingWhaleIdx)
                diveDepthRef(flyingWhaleIdx) = diveDepth_mean(itr_n)...
                    + diveDepth_std(itr_n)*randn(size(flyingWhaleIdx)); % add variation to dive depth,
                flyingWhaleIdx = find(diveDepthRef>=sd | diveDepthRef<1);
            end
            % Assign last n% to a descent phase
            % Choose a depth between start of clicking and destination depth
            % determine off-axis angle
            descentIdx = (floor((1-descentPerc(itr_n,1))*length(rho2))+1:length(rho2))';
            dFactor = rand(size(descentIdx));
            clickStartVec = clickStart_mean(itr_n,1) + clickStart_std(itr_n,1).*randn(size(descentIdx));
            
            % Find points that are above surface or below bottom and correct them
            flyingWhaleIdx2 = find(clickStartVec>=sd | clickStartVec<1);
            while ~isempty(flyingWhaleIdx2)
                clickStartVec(flyingWhaleIdx2) = clickStart_mean(itr_n,1)...
                    + clickStart_std(itr_n,1).*randn(size(flyingWhaleIdx2)); % add variation to dive depth,
                flyingWhaleIdx2 = find(clickStartVec>=sd | clickStartVec<1);
            end
            
            descentDelta = dFactor.* (diveDepthRef(descentIdx,:) - clickStartVec);
            diveDepthRef(descentIdx,1) = clickStartVec + descentDelta;
            
            %%%%% Beam Angle Computation %%%%%
            % Assign random beam orientation in horizontal (all orientations equally likely)
            randAngleVec = rand(size(rho2)).*359;
            theta2deg = theta2.*180./pi;
            partAngle = 180 + make180(thetaDeg);
            totalOffAxisHoriz = make180(randAngleVec - partAngle');
            onAxisHorz = abs(totalOffAxisHoriz) <= rotHorizDeg(itr_n,1) + ...
                rotHorizDeg_std(itr_n,1).*randn(1,size(totalOffAxisHoriz,2));
            
            % Compute vertical component of shift between animal and sensor (sd =
            % sensor depth)
            dZ = abs(sd - diveDepthRef);
            zAngle_180 = ceil(abs(atand(dZ./rho2)));
            zAngle_180(descentIdx,1) = ceil(abs(atand(dZ(descentIdx,:)./radRef(descentIdx,:))) - ...
            	 descAngle_mean(itr_n,1) + (descAngle_std(itr_n,1).*randn(size(descentIdx,1),1)));
            
            % If they're foraging, allow one amount of vertical rotation
            % if diving, different amount
            onAxisVert = abs(zAngle_180) <= rotVertDegForage(itr_n,1) + ...
                rotVertDegForageStd(itr_n,1)*randn(size(zAngle_180,1),1);
            
            onAxisVert(descentIdx,1)  = abs(zAngle_180(descentIdx,1) ) <= ...
            	(rotVertDegDive(itr_n,1) + rotVertDegDiveStd(itr_n,1)*randn(size(descentIdx,1),1));
            
            % Compute variation to add to source level
            SL_adj = SL_mean(itr_n,1) + (SL_std(itr_n,1)*randn(N,1));
            % directVec = directivity_mean(itr_n,1)+(directivity_std(itr_n,1)*randn(N,1));
            minAmp = minBeamAmp_mean(itr_n,1)+(minBeamAmp_std(itr_n,1)*randn(N,1));
            %%%%% Transmission Loss Loop %%%%%
            % initialize some variables
            RL = nan(size(thetaDeg));
            isheard = zeros(size(thetaDeg));
            distTL = nan(size(thetaDeg));
            
            for itr2 = 1:length(thetaDeg)
                % Compute location of this animal in the transmission loss matrix:
                % Find which row you want to look at:
                thisRd = rd_all{angleRef(itr2)};
                [~,thisDepthIdx] = min(abs(thisRd - round(diveDepthRef(itr2))));
                
                % Record the distance related portion of this transmission loss
                thisSortedTL = real(sortedTLVec{angleRef(itr2)});
                distTL(itr2) = thisSortedTL(thisDepthIdx,ceil(radRef(itr2)./rr_int));
                
                % Add up all the sources of TL
                if (onAxisVert(itr2)+onAxisHorz(itr2))==2
                    RL(itr2,1) = SL_adj(itr2) - distTL(itr2);
                else
                    RL(itr2,1) = SL_adj(itr2) - distTL(itr2) - minAmp(itr2);
                end
                % Is the total TL less than the maximum allowed?
                if thresh <= RL(itr2,1)
                    isheard(itr2,1) = 1; % detected it
                end
            end
            
            pDetTotal(itr_n,1) = sum(isheard)./length(isheard)';
            detVsLoc = [thetaDeg, rho2, isheard];
            totalSim = rho2';
            detSim = rho2(isheard==1)';
            
            % Compute detections in range bins, so you can make a histogram if desired
            % Makes more sense for click-based model
            % preallocate
            binTot = zeros(length(binVec)-1,1);
            binDet = zeros(length(binVec)-1,1);
            
            for itr2 = 1:length(binVec)-1
                binTot(itr2) = length(find(totalSim>binVec(itr2) & totalSim<binVec(itr2 +1)));
                binDet(itr2) = length(find(detSim>binVec(itr2) & detSim<binVec(itr2 +1)));
            end
            thisPercent = binDet./binTot;
            % save the bin counts to the overall set, so you can get means and variances per bin.
            binnedPercDet(itr_n,:) = thisPercent';
            binnedCountsDet(itr_n,:) = binDet';
            binnedCountsTot(itr_n,:) = binTot';
            
            
        end
        
        save(fullfile(saveDir,sprintf('%s_grpModOut_%dItrs_%s.mat',site,itr_n,species)),'-mat')
        % end
        
        % Histogram of detectability as a function of range
        spots = binVec(1:end-1)+(50);
        areas = ((binVec(2:end).^2)*pi)-((binVec(1:end-1).^2)*pi);
        means = nanmean(binnedPercDet)*100;
        means_keep = (means>0);
        spotsprune = spots(means_keep);
        means = means(means_keep);
        areas = areas(means_keep);
        errsTop = nanstd(binnedPercDet(:,means_keep)*100);
        errsBot = errsTop;
        toobig = (errsTop + means)>100;
        toosmall = (means - errsBot)<0;
        errsTop(toobig) = 100-means(toobig);
        errsBot(toosmall) = -(0-means(toosmall));
        
        figure(1); clf
        hb1 = bar(spotsprune,means,1);
        set(hb1,'EdgeColor','k','FaceColor','w')
        hold on
        ha = errorbar(spotsprune,means,errsBot,errsTop,'k');
        Xdata = get(ha,'Xdata');
        temp = 4:3:length(Xdata);
        temp(3:3:end) = [];
        % xleft and xright contain the indices of the left and right endpoints of the horizontal lines
        xleft = temp; xright = temp+1;
        Xdata(xleft) = Xdata(xleft) + 20;
        Xdata(xright) = Xdata(xright) - 20;
        set(ha,'Xdata',Xdata)
        plot(spotsprune,means,'-k','LineWidth',3)
        set(gca,'XTick',binVec,'FontSize',12)
        set(gca,'XTickLabel',binVec)
        xlabel(gca,'Horizontal Range (m)','FontSize',14)
        ylabel(gca, 'Probability of Detection (%)','FontSize',14)
        title({sprintf('Max Horiz. Range = %dm; mean P(det) = %1.2f%%; std = %1.2f%%', ...
            maxRange, nanmean(pDetTotal)*100, nanstd(pDetTotal)*100)},'FontSize',14)
        print(gcf,'-dpng','-r600',fullfile(saveDir,[site,'_',species,'_grpMod_pDet.png']))
        saveas(gca,fullfile(saveDir,[site,'_',species,'_grpMod_pDet.fig']))
        
        
        figure(2); clf
        binCountDetMean = nanmean(binnedCountsDet./...
            repmat(sum(binnedCountsDet,2),1,size(binnedCountsDet,2)));
        binCountDetStd = nanstd(binnedCountsDet./...
            repmat(sum(binnedCountsDet,2),1,size(binnedCountsDet,2)));
        binId = find(binCountDetMean-binCountDetStd<0);
        binCountStdBot = binCountDetStd;
        binCountStdBot(binId) = binCountDetMean(binId) ;
        hb2 = bar(spots/1000,binCountDetMean,1);
        hold on
        errorbar(spots/1000,binCountDetMean,binCountStdBot,binCountDetStd,'.k');
        set(hb2,'FaceColor',[1,1,1])
        set(gca,'XTick',binVec,'FontSize',12)
        set(gca,'XTickLabel',binVec)
        xlabel(gca,'Horizontal Range (km)','FontSize',14)
        ylabel(gca, '% of detections','FontSize',14)
        title('Kogia Group: Modeled Detection Ranges','FontSize',14)
        saveas(gca,fullfile(saveDir,[site,'_',species,'_grpMod_nDets.fig']))
        saveas(gca,fullfile(saveDir,[site,'_',species,'_grpMod_nDets.png']))
        
        % figure(3);clf
        % nDetsNormMean = mean(nDets(:,2:end))./mean(binnedCountsTot);
        % hb2 = bar(binVec(2:end)-50,nDetsNormMean,1);
        % set(hb2,'EdgeColor','k','FaceColor','w')
        % hold on
        % plot(binVec(2:end)-50,nDetsNormMean,'-k','LineWidth',3)
        % % errorbar(binVec(2:end)-50,nDetsNormMean,nDetsNormStdBot,nDetsNormStd,'*r')
        % set(gca,'XTick',binVec(1:5:end))
        % set(gca,'XTickLabel',binVec(1:5:end)/1000)
        % xlabel(gca,'Estimated Horizontal Range (km)','FontSize',14)
        % ylabel(gca, 'P(detection)','FontSize',14)
        % title({polarFile, 'Estimated P(detection) as a function of range based on recieved level'},'interpreter','none')
        % xlim([binVec(1),binVec(end)])
        % print(gcf,'-dpng','-r600',['E:\Data\John Reports\MC_kogia\',site,'_',species,'_grpMod_distanceEst.png'])
        % saveas(gca,[saveDir,site,'_',species,'_grpMod_distanceEst.fig'])
        %
        
        
        %
        % figure(4);clf
        % nDetsMean = mean(nDets(:,2:end),1);
        % nDetsStd = std(nDets(:,2:end),1);
        % binId3 = find(nDetsMean-nDetsStd<0);
        % nDetsStdBot = nDetsStd;
        % nDetsStdBot(binId3) = nDetsMean(binId3) ;
        %
        %
        % bar(binVec(2:end)-50,nDetsMean)
        % hold on
        %
        % errorbar(binVec(2:end)-50,nDetsMean,nDetsStdBot,nDetsStd,'*r')
        % set(gca,'XTick',binVec(1:5:end))
        % set(gca,'XTickLabel',binVec(1:5:end)/1000)
        % xlabel(gca,'Estimated Horizontal Range (km)','FontSize',14)
        % ylabel(gca, 'Counts','FontSize',14)
        % title({polarFile, 'Estimated Horizontal range based on recieved level'},'interpreter','none')
        % xlim([binVec(1),binVec(end)])
        % saveas(gca,[saveDir,site,'_',species,'_grpMod_distanceEst_counts.fig'])
        % saveas(gca,[saveDir,site,'_',species,'_grpMod_distanceEst_counts.png'])
    end
end