Major update on CCREClass and some bug fixes in the CRE class.

CCRE class is now "more logical" from a programing point of view.
Created a fucntion to create the bin-edges for the reference function.
"condensed" the calc function (removed redundancies)

CCREClass.m

@@ -12,10 +12,14 @@ classdef CCREClass<CREClass
     properties
         DataRef %This parameter will contain the reference set.
         nBinRef %number of "bins" to use for the reference distribution
+        EqualSizeBinFlagRef=false;  %set to true to get equally spaced bins, default : false, put an equal number of points per bin. 
+        EdgesRef; %edges for the bins of the reference function.  %todo, turn into dep. variable!
+    end
+    properties(SetAccess=private) %make read only
                 CCRE %parameter containing the conditional Cumulative Residual Entropy
-        EquidistantBinFlag=true; % set to true to get equidistant bins, set to false to get equi-content bins, 
+    end
-        % if set to false: each bin will contain (aproximately) an equal
+    properties(Constant)
-        % number of pixel values. 
+        UseHistProxyFlagRef=true; % this value is there to remind the user that we will have to bin the reference image. The actual value is not used in the code
     end
     
     methods %constructor
@@ -41,6 +45,16 @@ classdef CCREClass<CREClass
         end
     end
     methods %user callable functions
+        function CalcEdgesRef(obj)
+            % copy the info for "ref" to a CRE object
+            Stmp=CREClass('ref',obj.DataRef...
+                ,'nbin',obj.nBinRef...
+                ,'EqualSizeBinFlag',obj.EqualSizeBinFlagRef...
+                ,'UseHistProxyFlag',obj.UseHistProxyFlagRef...
+                );
+            Stmp.CalcEdges;
+            obj.EdgesRef=Stmp.Edges;
+        end
         %% CCRE A|B
         %
         % $$CCRE = \varepsilon (A) - E[\varepsilon(A|B)]$$
@@ -61,21 +75,18 @@ classdef CCREClass<CREClass
             if isempty(obj.DataRef);return;end %if no refdata given, return;
             %calculate histogram for RefData
             EVarEpsDataGivenRef=0; %clear previous calculations
-            if obj.EquidistantBinFlag
+            %get the edge definition. 
-                [pDataRef,E]=histcounts(obj.DataRef(:),obj.nBinRef,'Normalization','probability');
+            % a bit of a workaround, I don't plan to reprogram the function
-                [JointHist,E1,E2]=histcounts2(obj.Data(:),obj.DataRef(:),[obj.nBin,obj.nBinRef],'Normalization','probability'); %obtain bin edges for Target
+            if isempty(obj.EdgesRef)
-            else
+                obj.CalcEdgesRef;
-                prc=linspace(0,100,obj.nBinRef+1); % define the edges of the bins as a precentile
-                Edges=prctile(obj.DataRef(:),prc);
-                [pDataRef,E]=histcounts(obj.DataRef(:),Edges,'Normalization','probability');
-                [JointHist,E1,E2]=histcounts2(obj.Data(:),obj.DataRef(:),[obj.nBin,obj.nBinRef],'Normalization','probability'); %obtain bin edges for Target
             end
+            [pDataRef,E]=histcounts(obj.DataRef(:),obj.EdgesRef,'Normalization','probability');
             for k=1:numel(pDataRef) %loop over the bins in pDataRef NB: must be equal to the bins in JointHist for the Ref
                 if pDataRef(k)==0;continue;end% if no data at this line go to the next
                 ind=E(k)<=obj.DataRef(:)&obj.DataRef(:)<E(k+1); %get an index to all points in Ref data inside the current bin
-                CRE_DataGivenDataRef=CREClass;
+                CRE_DataGivenDataRef=CREClass('data',obj.Data(ind)...
-                CRE_DataGivenDataRef.Data=obj.Data(ind);
+                    ,'nbin',obj.Edges...
-                CRE_DataGivenDataRef.nBin=E1; %force the same bin definition for all itterations. This value is ignored if UseHistProxy is set to false
+                    );
                 CRE_DataGivenDataRef.Calc;
                 EVarEpsDataGivenRef=EVarEpsDataGivenRef+pDataRef(k)*CRE_DataGivenDataRef.CRE; % I really have to check this line on merit. 
             end

CREClass.m

@@ -3,7 +3,7 @@ classdef CREClass < handle
         Data; % the  data set to get the CRE off
         nBin; %number of bins 
         Monitor=true; %set to true to get intermediat results/figures
-        UseHistProxy=true; % set this to get a proxy of the survival function using the hist function
+        UseHistProxyFlag=true; % set this to get a proxy of the survival function using the hist function
         EqualSizeBinFlag=true; % set to true to get equal sized bins, set to false to get equi-content bins, 
         % if set to false: each bin will contain (aproximately) an equal
         % number of pixel values. 
@@ -12,11 +12,15 @@ classdef CREClass < handle
         CRE
         RB
         P_RB
+    end
+    
+    properties(Dependent)
                 Edges
     end
     
     properties (Hidden)
         WaitBarHandle;
+        HiddenEdges;
     end
     properties (Constant,Hidden)
         alpha=1; %needed for the scaled cummulative redual entropy as defined in zografoss. Will be implemented in the future. (I hope)
@@ -38,6 +42,10 @@ classdef CREClass < handle
                             obj.Data=varargin{k};
                         case {'nbin','bin'}
                             obj.nBin=varargin{k};
+                        case {'equalsizebinflag'}
+                            obj.EqualSizeBinFlag=varargin{k};
+                        case {'usehistproxyflag'}
+                            obj.UseHistProxyFlag=varargin{k};
                         otherwise
                             error('invalid variable name');
                     end
@@ -52,8 +60,17 @@ classdef CREClass < handle
             end
         end
     end
-    %% set functions
+    %% set and get functions
     methods
+        function Edges=get.Edges(obj)
+            Edges=obj.HiddenEdges;
+        end
+        function set.Edges(obj,val)
+            if ~isnumeric(val);error('provide a array with edge values');end
+            if numel(val)<3;error('too few edges provided. Need at least 3 values, i.e. 2 bins defined');end
+            obj.HiddenEdges=transpose(val(:));
+            obj.nBin=numel(val)-1;
+        end
         function set.Data(obj,val)
             if isnumeric(val)
                 obj.Data=val;
@@ -76,7 +93,7 @@ classdef CREClass < handle
                 warning('nothing to do');
                 return;
              end
-             if obj.UseHistProxy
+             if obj.UseHistProxyFlag
                 if isempty(obj.nBin)
                     obj.nBin=floor(numel(obj.Data)./10);
                 end
@@ -87,7 +104,7 @@ classdef CREClass < handle
                         obj.Edges=prctile(obj.Data(:),prc); %find the edges of the bins
                 end
              else
-                 if obj.UseHistProxy ; error('You cannot use EqualSizeBinFlag=true and UseHistProxy=false as a combination');end
+                 if obj.UseHistProxyFlag ; error('You cannot use EqualSizeBinFlag=true and UseHistProxyFlag=false as a combination');end
                  q=transpose(unique(sort(obj.Data(:))));
                  dq=diff(q);
                  obj.Edges=cat(2,q(1)-eps,q(1:end-1)+dq/2,q(end)+eps); %the eps is a "shortcut" to avoid value=edge cases. Note that eps will increase the size of the bin by a very small amount. Although this is a systemetic bias, its effect will be neglectable
@@ -118,7 +135,7 @@ classdef CREClass < handle
             dl(FC>0.5)=0; %set the weight for all in the histogram to the left side (i.e. <p50) to 0
             CRE_Med_Inf=-transpose(dl)*(FC.*LogFC)./(transpose(dl)*FC); %CRE zografos 
             ind=find(exp(-CRE_Med_Inf)>=FC,1,'first');
-            if obj.UseHistProxy
+            if obj.UseHistProxyFlag
                 obj.P_RB=FC(ind);
                 if ind==numel(bin)%if RB is (beyond) the last bin,
                     warning(sprintf('RB is hitting the upper bin value\n Maybe increase number of bins')) %#ok<*SPWRN>

TestCREClass.m

@@ -6,30 +6,30 @@ A=randn(100);
 %% create instance of CRE class
 S1=CREClass;
 %% Calculate by setting labda equal to each of the (unique) values in the data. 
-S1.UseHistProxy=false;
+S1.UseHistProxyFlag=false;
 S1.EqualSizeBinFlag=false;
 S1.Data=A;% set gaussian data
 S1.Calc;
-
+disp(S1)
-% use histogram aproximation with unequal bin sizes in the histogram
+%% use histogram aproximation with unequal bin sizes in the histogram
 S2=CREClass;
-S2.UseHistProxy=true;
+S2.UseHistProxyFlag=true;
 S2.EqualSizeBinFlag=false;
 S2.Data=A;% set gaussian data
 S2.nBin=numel(A);      % set nbin to npoints; This is fine as we use the cummulative distribution and the formulas as defined in Zografos
 S2.Calc
-
+disp(S2)
-% use histogram aproximation with equal bin sizes in the histogram
+%% use histogram aproximation with equal bin sizes in the histogram
 S3=CREClass;
-S3.UseHistProxy=true;
+S3.UseHistProxyFlag=true;
 S3.EqualSizeBinFlag=false;
 S3.Data=A;% set gaussian data
 S3.nBin=numel(A);      % set nbin to npoints; This is fine as we use the cummulative distribution and the formulas as defined in Zografos
 S3.Calc
-
+disp(S3)
-%Check recalculation of edges
+%% Check recalculation of edges
 S4=CREClass;
-S4.UseHistProxy=true;
+S4.UseHistProxyFlag=true;
 S4.EqualSizeBinFlag=true;
 S4.Data=A;% set gaussian data
 S4.nBin=numel(A);      % set nbin to npoints; This is fine as we use the cummulative distribution and the formulas as defined in Zografos
@@ -40,16 +40,16 @@ S4.CalcEdges ; % force a recalculation of the edges
 S4.Calc;
 disp(S4)
 
-return;
+
 %% Show effect of scaling or ofsett the data
 % Scale
 S.Data=100*A;
-S.UseHistProxy=true;
+S.UseHistProxyFlag=true;
 S.Calc;
 disp(S);
 % Offset
 S.Data=A+10;
-S.UseHistProxy=true;
+S.UseHistProxyFlag=true;
 S.Calc;
 disp(S);
 %% 
@@ -66,7 +66,7 @@ out(numel(w),1)=struct('w',[],'CRE',[],'CCRE',[],'R',[]);
 for k=1:numel(w)
     out(k).w=w(k);
     CS=CCREClass;
-    CS.EquidistantBinFlag=false;
+    CS.EqualSizeBinFlag=false;
     CS.Data=1*(1-abs(w(k)))*A+w(k)*B;
     CS.nBin=50;
     CS.DataRef=B;
@@ -86,4 +86,4 @@ disp(CS);
     plot([out(:).w],[out(:).CCRE],'o');
     plot([out(:).w],[out(:).CRE],'x');
     figure(3);clf
-    plot([out(:).w],[out(:).R],'d')
+    plot([out(:).w],[out(:).R],'d');