replaces randsample with a custom-built randSample function and removes range

examples/normalReflectivity/customXY/customXYDSPCScript.m

@@ -218,7 +218,7 @@
 
     % We don't need to calculate all of it, just take a random 1000 points
     % from the chain. Make a set of indices...
-    samples = randsample(nSamples,1000);
+    samples = randSample(nSamples,1000);
 
     % Make some empty arrays to store our data....
     vfSi = []; vfOxide = []; vfHeadL = []; vfTails = []; vfHeadR = []; vfWat = [];
@@ -281,4 +281,4 @@
 shade(z,ciHeadR(1,:),z,ciHeadR(2,:),'FillColor',cols(5,:),'FillType',fillType,'FillAlpha',fillAlpha);
 title('Volume Fractions');
 %% 
-% .. and we are done.
+% .. and we are done.

minimisers/DREAM/functions/calcProposal.m

@@ -5,7 +5,7 @@
 % % % eps = DREAMPar.zeta * randn(DREAMPar.nChains,DREAMPar.nParams);
 % % % 
 % % % % Determine which sequences to evolve with what DE strategy
-% % % DE_pairs = randsample( [1:DREAMPar.delta ] , DREAMPar.nChains , true , [ 1/DREAMPar.delta*ones(1,DREAMPar.delta) ])';
+% % % DE_pairs = randSample( [1:DREAMPar.delta ] , DREAMPar.nChains , [ 1/DREAMPar.delta*ones(1,DREAMPar.delta) ])';
 % % % 
 % % % % Generate series of permutations of chains
 % % % [dummy,tt] = sort(rand(DREAMPar.nChains-1,DREAMPar.nChains));
@@ -100,7 +100,7 @@
 eps = DREAMPar.zeta * randn(DREAMPar.nChains,DREAMPar.nParams);
 
 % Determine how many chain pairs to use for each individual chain
-DE_pairs = randsample( [1:DREAMPar.delta ] , DREAMPar.nChains , true , [ 1/DREAMPar.delta*ones(1,DREAMPar.delta) ])';
+DE_pairs = randSample( [1:DREAMPar.delta ] , DREAMPar.nChains , [ 1/DREAMPar.delta*ones(1,DREAMPar.delta) ])';
 
 % Generate uniform random numbers for each chain to determine which dimension to update
 rnd_cr = rand(DREAMPar.nChains,DREAMPar.nParams);
@@ -116,7 +116,7 @@
 dx = zeros(DREAMPar.nChains,DREAMPar.nParams);               
 
 % Determine when jumprate is 1
-gamma = randsample([0 1],DREAMPar.nChains,true,[ 1 - DREAMPar.pUnitGamma DREAMPar.pUnitGamma ]);   
+gamma = randSample([0 1],DREAMPar.nChains,[ 1 - DREAMPar.pUnitGamma DREAMPar.pUnitGamma ]);   
 
 % Create N proposals
 for i = 1:DREAMPar.nChains               
@@ -153,4 +153,4 @@
 % If specified do boundary handling ( "Bound","Reflect","Fold")
 if isfield(paramInfo,'boundhandling')
     [x_new] = boundaryHandling(x_new,paramInfo);
-end
+end

minimisers/DREAM/functions/drawCR.m

@@ -35,6 +35,6 @@
         CR = reshape(cCR,DREAMPar.nChains,[]);
     else
         % If crossover probabilities are not updated  
-        CR = reshape(randsample([1:DREAMPar.nCR]/DREAMPar.nCR,DREAMPar.steps*DREAMPar.nChains,true,pCR),DREAMPar.nChains,DREAMPar.steps);
+        CR = reshape(randSample([1:DREAMPar.nCR]/DREAMPar.nCR,DREAMPar.steps*DREAMPar.nChains,pCR),DREAMPar.nChains,DREAMPar.steps);
     end    
 end

utilities/misc/randSample.m

@@ -0,0 +1,57 @@
+function outputSample = randSample(population, numItems, weights)
+  % Take a random sample of values.
+  % GPL-licensed replacement for `randsample` from Statistics and Machine
+  % Learning Toolbox.
+  % 
+  % Note that unweighted sampling is *without* replacement, and weighted
+  % sampling is *with* replacement!
+  % 
+  % Parameters
+  % ----------
+  % population : vector or int
+  %   if a vector, sample k values from the vector.
+  %   if an int, sample k values from 1:n.
+  % numItems : int
+  %   the number of items to sample.
+  % weights : vector
+  %   a weight vector, where the i'th index of ``weights`` is the 
+  %   weight for the i'th member of the population.
+
+  % if the user is trying to get multiple items from a single-integer array,
+  % assume that they want to sample between 1 and n
+  if isscalar(population) && numItems > 1
+    population = 1:population;
+  end
+
+  nargs = nargin;
+
+  if nargs > 3
+    error("Too many arguments.")
+  end
+
+  if nargs == 3
+
+    if length(weights) ~= length(population)
+      error("Weights and population must be the same length.")
+    end
+
+    % we generate weighted random integers by creating bins from our weights
+    % and discretizing random numbers in [0, 1] to those bins
+    weights = normalize(weights, 'norm', 1);
+    bins = [0, cumsum(weights)];
+    bins(end) = 1;  % ensure final bin is 1 as normalize is not always exact
+
+    randomIndices = discretize(rand(1, numItems), bins);
+
+    outputSample = population(randomIndices);
+
+  % we can generate unweighted numbers far more efficiently
+  % just randomise the array and return the first numItems items
+  else
+    if numItems > length(population)
+      error("numItems is larger than the number of items in the population.")
+    end
+      population = population(randperm(length(population)));
+      outputSample = population(1:numItems);
+  end
+end

utilities/plotting/cornerPlot.m

@@ -308,7 +308,7 @@ function smoothhist2DLocal(X,lambda,nbins,outliercutoff,plottype)
 %         plot(X(outliers,1),X(outliers,2),'.','MarkerEdgeColor',[.8 .8 .8]);
 %     end
 % %     % plot a subsample of the data
-% %     Xsample = X(randsample(n,n/10),:);
+% %     Xsample = X(randSample(n,n/10),:);
 % %     plot(Xsample(:,1),Xsample(:,2),'bo');
 %     hold off
 % end

utilities/plotting/histp.m

@@ -35,7 +35,7 @@
 if nargin<2|isempty(nbin)
 % nbin=ceil(log2(n)+1);                        % Sturges' formula
 % nbin=ceil(range(y)/(3.5*sy*n^(-1/3)));       % Scott
-  nbin=ceil(range(y)/(2*iqrange(y)*n^(-1/3))); % Freedman & Diaconis
+  nbin=ceil((max(y) - min(y))/(2*iqrange(y)*n^(-1/3))); % Freedman & Diaconis
 
   % bar(...,'w') does not work with nbin>150
   nbin=min(50,nbin);