function [y,ax,bx] = scaleCols(x, x_scale, optional_bx, optional_action)
%  Scales a column matrix of observations so that each row will occupy
%  the given interval [ x_scale(1) x_scale(2) ]. Alternatively re-uses 
%  or inverts a former scaling.
%
%  If I'm not mistaken, this originates from a course excercise set by
%  Erkki Heikkola (for the course "Application programming in Matlab")
%
%  There are three ways to call this:
%
%     [y,ax,bx] = scaledata(x, x_scale) 
%         analyzes the data in x to find the suitable transformation and
%         performs the scaling
%
%
%     y = scaledata(x, ax, bx) 
%         re-uses a formerly found transformation. In the neural computing
%         territory, you would ordinarily create the scaling transformation
%         according to your known training data prior to training. When
%         applying the network to previously unseen data, the scaling must
%         be congruent with the one used in training; otherwise there would
%         be systematically wrong results due to mismatched input
%         magnitudes.
%
%     x = scaledata(y, ax, bx, 'inverse')
%         re-uses a transformation by providing its inverse.
%
%  
%  Attention: 
%
%    - each row of the matrix x corresponds to
%      a specific variable and each column to a specific observation.
%      The purpose of this function is to match the range of variables to
%      the range most suitable for, for example, feed-forward neural 
%      network.
%
%    - if some variable has a near-constant value according to the 
%      observations, its value will be scaled to the middle of the given 
%      interval
%
%    - to use a previously discovered scaling for a new data set z, you can 
%      call this routine with the ax and bx parameters, or write:
%
%         scaled_z = diag(ax)*z + repmat(bx,1,size(z,2));
%
% 
%  Input
%    x         matrix to be scaled
%    x_scale   interval as a two-vector.
%
%  Output
%    y         scaled matrix
%    ax        information about the scaling;
%    bx        required for inverse transformation
%

% Action 1, "re-scaling":
if ((nargin == 3) && (nargout == 1))
    ax = x_scale; % a naming oddity due to "DIY-overloading".
    bx = optional_bx;
    y = diag(ax)*x + repmat(bx,1,size(x,2));
    return
end

% Action 2, "inverse-scaling":
if ((nargin == 4) && (nargout == 1) && strcmp(optional_action,'inverse'))
    ax = x_scale; % a naming oddity due to "DIY-overloading".
    bx = optional_bx;
    y = (x - repmat(bx,1,size(x,2))) ./ repmat(ax,1,size(x,2));
    return
end


% Action 3, "scaling based on analyzing the input":
if (x_scale(2) - x_scale(1)) <= 1.e-13
   error('scaledata: Small or inverted scaling intervals not allowed.');
end


a   = min(x,[],2);
b   = max(x,[],2);
bma = b - a;

ind = find(bma <= 1.e-13);
bma(ind) = x_scale(2) - x_scale(1);

ax  = (x_scale(2) - x_scale(1)) ./ bma;
bx  = x_scale(1) - ax .* a;
bx(ind) = (x_scale(1)+x_scale(2))/2 - a(ind);

y = diag(ax)*x + repmat(bx,1,size(x,2));
