% runica() - Perform Independent Component Analysis (ICA) decomposition
%            of input data using the logistic infomax ICA algorithm of
%            Bell & Sejnowski (1995) with the natural gradient feature
%            of Amari, Cichocki & Yang, or optionally the extended-ICA
%            algorithm of Lee, Girolami & Sejnowski, with optional PCA
%            dimension reduction. Annealing based on weight changes is
%            used to automate the separation process.
% Usage:
%         >> [weights,sphere] = runica(data); % train using defaults
%    else
%         >> [weights,sphere,compvars,bias,signs,lrates,activations] ...
%                             = runica(data,'Key1',Value1',...);
% Input:
%    data     = input data (chans,frames*epochs).
%               Note that if data consists of multiple discontinuous epochs,
%               each epoch should be separately baseline-zero'd using
%                  >> data = rmbase(data,frames,basevector);
%
% Optional keywords [argument]:
% 'extended'  = [N] perform tanh() "extended-ICA" with sign estimation
%               N training blocks. If N > 0, automatically estimate the
%               number of sub-Gaussian sources. If N < 0, fix number of
%               sub-Gaussian comps to -N [faster than N>0] (default|0 -> off)
% 'pca'       = [N] decompose a principal component     (default -> 0=off)
%               subspace of the data. Value is the number of PCs to retain.
% 'sphering'  = ['on'/'off'] flag sphering of data      (default -> 'on')
% 'weights'   = [W] initial weight matrix               (default -> eye())
%                            (Note: if 'sphering' 'off', default -> spher())
% 'lrate'     = [rate] initial ICA learning rate (<< 1) (default -> heuristic)
% 'block'     = [N] ICA block size (<< datalength)      (default -> heuristic)
% 'anneal'    = annealing constant (0,1] (defaults -> 0.90, or 0.98, extended)
%                         controls speed of convergence
% 'annealdeg' = [N] degrees weight change for annealing (default -> 70)
% 'stop'      = [f] stop training when weight-change < this (default -> 1e-6
%               if less than 33 channel and 1E-7 otherwise)
% 'maxsteps'  = [N] max number of ICA training steps    (default -> 512)
% 'bias'      = ['on'/'off'] perform bias adjustment    (default -> 'on')
% 'momentum'  = [0<f<1] training momentum               (default -> 0)
% 'specgram'  = [srate loHz hiHz frames winframes] decompose a complex time/frequency
%               transform of the data - though not optimally. (Note: winframes must
%               divide frames) (defaults [srate 0 srate/2 size(data,2) size(data,2)])
% 'posact'    = make all component activations net-positive(default 'off'}
%               Requires time and memory; posact() may be applied separately.
% 'verbose'   = give ascii messages ('on'/'off')        (default -> 'on')
% 'logfile'   = [filename] save all message in a log file in addition to showing them
%               on screen (default -> none)
% 'interrupt' = ['on'|'off'] draw interrupt figure. Default is off.
% 'rndreset'  = ['on'|'off'] reset the random seed. Default is off (although it
%               used to be on prior to 2015. This means that ICA will always return.
%               the same decomposition unless this option is set to 'on'.
%
% Outputs:    [Note: RO means output in reverse order of projected mean variance
%                    unless starting weight matrix passed ('weights' above)]
% weights     = ICA weight matrix (comps,chans)      [RO]
% sphere      = data sphering matrix (chans,chans) = spher(data)
%               Note that unmixing_matrix = weights*sphere {if sphering off -> eye(chans)}
% compvars    = back-projected component variances   [RO]
% bias        = vector of final (ncomps) online bias [RO]    (default = zeros())
% signs       = extended-ICA signs for components    [RO]    (default = ones())
%                   [ -1 = sub-Gaussian; 1 = super-Gaussian]
% lrates      = vector of learning rates used at each training step [RO]
% activations = activation time courses of the output components (ncomps,frames*epochs)
%
% Authors: Scott Makeig with contributions from Tony Bell, Te-Won Lee,
% Tzyy-Ping Jung, Sigurd Enghoff, Michael Zibulevsky, Delorme Arnaud,
% CNL/The Salk Institute, La Jolla, 1996-

% Uses: posact()

% 'ncomps'    = [N] number of ICA components to compute (default -> chans or 'pca' arg)
%               using rectangular ICA decomposition. This parameter may return
%               strange results. This is because the weight matrix is rectangular
%               instead of being square. Do not use except to try to fix the problem.

% Reference (please cite):
%
% Makeig, S., Bell, A.J., Jung, T-P and Sejnowski, T.J.,
% "Independent component analysis of electroencephalographic data,"
% In: D. Touretzky, M. Mozer and M. Hasselmo (Eds). Advances in Neural
% Information Processing Systems 8:145-151, MIT Press, Cambridge, MA (1996).
%
% Toolbox Citation:
%
% Makeig, Scott et al. "EEGLAB: ICA Toolbox for Psychophysiological Research".
% WWW Site, Swartz Center for Computational Neuroscience, Institute of Neural
% Computation, University of San Diego California
% <www.sccn.ucsd.edu/eeglab/>, 2000. [World Wide Web Publication].
%
% For more information:
% http://www.sccn.ucsd.edu/eeglab/icafaq.html - FAQ on ICA/EEG
% http://www.sccn.ucsd.edu/eeglab/icabib.html - mss. on ICA & biosignals
% http://www.cnl.salk.edu/~tony/ica.html - math. mss. on ICA

% Copyright (C) 1996 Scott Makeig et al, SCCN/INC/UCSD, scott@sccn.ucsd.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

%%%%%%%%%%%%%%%%%%%%%%%%%%% Edit history %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  runica()  - by Scott Makeig with contributions from Tony Bell, Te-Won Lee
%              Tzyy-Ping Jung, Sigurd Enghoff, Michael Zibulevsky et al.
%                            CNL / Salk Institute 1996-00
%  04-30-96 built from icatest.m and ~jung/.../wtwpwica.m -sm
%  07-28-97 new runica(), adds bias (default on), momentum (default off),
%           extended-ICA (Lee & Sejnowski, 1997), cumulative angledelta
%           (until lrate drops), keywords, signcount for speeding extended-ICA
%  10-07-97 put acos() outside verbose loop; verbose 'off' wasn't stopping -sm
%  11-11-97 adjusted help msg -sm
%  11-30-97 return eye(chans) if sphering 'off' or 'none' (undocumented option) -sm
%  02-27-98 use pinv() instead of inv() to rank order comps if ncomps < chans -sm
%  04-28-98 added 'posact' and 'pca' flags  -sm
%  07-16-98 reduced length of randperm() for kurtosis subset calc. -se & sm
%  07-19-98 fixed typo in weights def. above -tl & sm
%  12-21-99 added 'specgram' option suggested by Michael Zibulevsky, UNM -sm
%  12-22-99 fixed rand() sizing inefficiency on suggestion of Mike Spratling, UK -sm
%  01-11-00 fixed rand() sizing bug on suggestion of Jack Foucher, Strasbourg -sm
%  12-18-00 test for existence of Sig Proc Tlbx function 'specgram'; improve
%           'specgram' option arguments -sm
%  01-25-02 reformated help & license -ad
%  01-25-02 lowered default lrate and block -ad
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [weights,sphere,meanvar,bias,signs,lrates,data,y] = runica(data,varargin) % NB: Now optionally returns activations as variable 'data' -sm 7/05

if nargin < 1
  help runica
  return
end

[chans frames] = size(data); % determine the data size
urchans = chans;  % remember original data channels
datalength = frames;
if chans<2
   fprintf('\nrunica() - data size (%d,%d) too small.\n\n', chans,frames);
   return
end
%
%%%%%%%%%%%%%%%%%%%%%% Declare defaults used below %%%%%%%%%%%%%%%%%%%%%%%%
%
MAX_WEIGHT           = 1e8;       % guess that weights larger than this have blown up
DEFAULT_STOP         = 0.000001;  % stop training if weight changes below this
DEFAULT_ANNEALDEG    = 60;        % when angle change reaches this value,
DEFAULT_ANNEALSTEP   = 0.90;      %     anneal by multiplying lrate by this
DEFAULT_EXTANNEAL    = 0.98;      %     or this if extended-ICA
DEFAULT_MAXSTEPS     = 512;       % ]top training after this many steps
DEFAULT_MOMENTUM     = 0.0;       % default momentum weight

DEFAULT_BLOWUP       = 1000000000.0;   % = learning rate has 'blown up'
DEFAULT_BLOWUP_FAC   = 0.8;       % when lrate 'blows up,' anneal by this fac
DEFAULT_RESTART_FAC  = 0.9;       % if weights blowup, restart with lrate
                                  % lower by this factor
MIN_LRATE            = 0.000001;  % if weight blowups make lrate < this, quit
MAX_LRATE            = 0.1;       % guard against uselessly high learning rate
DEFAULT_LRATE        = 0.00065/log(chans);
                                  % heuristic default - may need adjustment
                                  %   for large or tiny data sets!
% DEFAULT_BLOCK        = floor(sqrt(frames/4));  % heuristic default
DEFAULT_BLOCK          = ceil(min(5*log(frames),0.3*frames)); % heuristic
                                  % - may need adjustment!
% Extended-ICA option:
DEFAULT_EXTENDED     = 0;         % default off
DEFAULT_EXTBLOCKS    = 1;         % number of blocks per kurtosis calculation
DEFAULT_NSUB         = 1;         % initial default number of assumed sub-Gaussians
                                  % for extended-ICA
DEFAULT_EXTMOMENTUM  = 0.5;       % momentum term for computing extended-ICA kurtosis
MAX_KURTSIZE         = 6000;      % max points to use in kurtosis calculation
MIN_KURTSIZE         = 2000;      % minimum good kurtosis size (flag warning)
SIGNCOUNT_THRESHOLD  = 25;        % raise extblocks when sign vector unchanged
                                  % after this many steps
SIGNCOUNT_STEP       = 2;         % extblocks increment factor

DEFAULT_SPHEREFLAG   = 'on';      % use the sphere matrix as the default
                                  %   starting weight matrix
DEFAULT_INTERRUPT    = 'off';     % figure interruption
DEFAULT_PCAFLAG      = 'off';     % don't use PCA reduction
DEFAULT_POSACTFLAG   = 'off';     % don't use posact(), to save space -sm 7/05
DEFAULT_VERBOSE      = 1;         % write ascii info to calling screen
DEFAULT_BIASFLAG     = 1;         % default to using bias in the ICA update rule
DEFAULT_RESETRANDOMSEED = true;   % default to reset the random number generator to a 'random state'

%
%%%%%%%%%%%%%%%%%%%%%%% Set up keyword default values %%%%%%%%%%%%%%%%%%%%%%%%%
%
if nargout < 2,
    fprintf('runica() - needs at least two output arguments.\n');
    return
end
epochs = 1;							 % do not care how many epochs in data

pcaflag    = DEFAULT_PCAFLAG;
sphering   = DEFAULT_SPHEREFLAG;     % default flags
posactflag = DEFAULT_POSACTFLAG;
verbose    = DEFAULT_VERBOSE;
logfile    = [];

block      = DEFAULT_BLOCK;          % heuristic default - may need adjustment!
lrate      = DEFAULT_LRATE;
annealdeg  = DEFAULT_ANNEALDEG;
annealstep = 0;                      % defaults declared below
nochange   = NaN;
momentum   = DEFAULT_MOMENTUM;
maxsteps   = DEFAULT_MAXSTEPS;

weights    = 0;                      % defaults defined below
ncomps     = chans;
biasflag   = DEFAULT_BIASFLAG;

interrupt  = DEFAULT_INTERRUPT;
extended   = DEFAULT_EXTENDED;
extblocks  = DEFAULT_EXTBLOCKS;
kurtsize   = MAX_KURTSIZE;
signsbias  = 0.02;                   % bias towards super-Gaussian components
extmomentum= DEFAULT_EXTMOMENTUM;    % exp. average the kurtosis estimates
nsub       = DEFAULT_NSUB;
wts_blowup = 0;                      % flag =1 when weights too large
wts_passed = 0;                      % flag weights passed as argument
reset_randomseed = DEFAULT_RESETRANDOMSEED;

%
%%%%%%%%%% Collect keywords and values from argument list %%%%%%%%%%%%%%%
%
   if (nargin> 1 && rem(nargin,2) == 0)
      fprintf('runica(): Even number of input arguments???')
      return
   end
   for i = 1:2:length(varargin) % for each Keyword
      Keyword = varargin{i};
      Value = varargin{i+1};
      if ~ischar(Keyword)
         fprintf('runica(): keywords must be strings')
         return
      end
      Keyword = lower(Keyword); % convert upper or mixed case to lower

      if strcmp(Keyword,'weights') || strcmp(Keyword,'weight')
         if ischar(Value)
            fprintf(...
      'runica(): weights value must be a weight matrix or sphere')
            return
         else
           weights = Value;
           wts_passed =1;
         end
      elseif strcmp(Keyword,'ncomps')
         if ischar(Value)
            fprintf('runica(): ncomps value must be an integer')
            return
         end
         if ncomps < urchans && ncomps ~= Value
            fprintf('runica(): Use either PCA or ICA dimension reduction');
            return
         end
         fprintf('*****************************************************************************************');
         fprintf('************** WARNING: NCOMPS OPTION OFTEN DOES NOT RETURN ACCURATE RESULTS ************');
         fprintf('************** WARNING: IF YOU FIND THE PROBLEM, PLEASE LET US KNOW          ************');
         fprintf('*****************************************************************************************');
         ncomps = Value;
         if ~ncomps,
            ncomps = chans;
         end
      elseif strcmp(Keyword,'pca')
         if ncomps < urchans && ncomps ~= Value
            fprintf('runica(): Use either PCA or ICA dimension reduction');
            return
         end
         if ischar(Value)
            fprintf(...
'runica(): pca value should be the number of principal components to retain')
            return
         end
         pcaflag = 'on';
         ncomps = Value;
         if ncomps > chans || ncomps < -1
            fprintf('runica(): pca value must be in range [%d,%d]\n',-chans+1, chans)
            return
         end
         if ncomps < 0
             ncomps = size(data,1)+ncomps;
         end
         chans = ncomps;
       elseif strcmp(Keyword,'interupt') || strcmp(Keyword,'interrupt')
         if ~ischar(Value)
           fprintf('runica(): interrupt value must be on or off')
           return
         else
           Value = lower(Value);
           if ~strcmp(Value,'on') && ~strcmp(Value,'off')
             fprintf('runica(): interrupt value must be on or off')
             return
           end
           interrupt = Value;
         end
     elseif strcmp(Keyword,'posact')
         if ~ischar(Value)
           fprintf('runica(): posact value must be on or off')
           return
         else
           Value = lower(Value);
           if ~strcmp(Value,'on') && ~strcmp(Value,'off'),
             fprintf('runica(): posact value must be on or off')
             return
           end
           posactflag = Value;
         end
      elseif strcmp(Keyword,'lrate')
         if ischar(Value)
            fprintf('runica(): lrate value must be a number')
            return
         end
         lrate = Value;
         if lrate>MAX_LRATE || lrate <0
           fprintf('runica(): lrate value is out of bounds');
           return
         end
         if ~lrate
            lrate = DEFAULT_LRATE;
         end
      elseif strcmp(Keyword,'block') || strcmp(Keyword,'blocksize')
         if ischar(Value)
            fprintf('runica(): block size value must be a number')
            return
         end
         block = floor(Value);
         if ~block
           block = DEFAULT_BLOCK;
         end
      elseif strcmp(Keyword,'stop') || strcmp(Keyword,'nochange') ...
                    || strcmp(Keyword,'stopping')
         if ischar(Value)
            fprintf('runica(): stop wchange value must be a number')
            return
         end
         nochange = Value;
      elseif strcmp(Keyword,'logfile')
         if ~ischar(Value)
            fprintf('runica(): logfile value must be a string')
            return
         end
         logfile = Value;
      elseif strcmp(Keyword,'maxsteps') || strcmp(Keyword,'steps')
         if ischar(Value)
            fprintf('runica(): maxsteps value must be an integer')
            return
         end
         maxsteps = Value;
         if ~maxsteps
            maxsteps   = DEFAULT_MAXSTEPS;
         end
         if maxsteps < 0
            fprintf('runica(): maxsteps value (%d) must be a positive integer',maxsteps)
            return
         end
      elseif strcmp(Keyword,'anneal') || strcmp(Keyword,'annealstep')
         if ischar(Value)
            fprintf('runica(): anneal step value (%2.4f) must be a number (0,1)',Value)
            return
         end
         annealstep = Value;
         if annealstep <=0 || annealstep > 1
            fprintf('runica(): anneal step value (%2.4f) must be (0,1]',annealstep)
            return
         end
      elseif strcmp(Keyword,'annealdeg') || strcmp(Keyword,'degrees')
         if ischar(Value)
            fprintf('runica(): annealdeg value must be a number')
            return
         end
         annealdeg = Value;
         if ~annealdeg
             annealdeg = DEFAULT_ANNEALDEG;
         elseif annealdeg > 180 || annealdeg < 0
          fprintf('runica(): annealdeg (%3.1f) is out of bounds [0,180]',...
                annealdeg);
          return

         end
      elseif strcmp(Keyword,'momentum')
         if ischar(Value)
            fprintf('runica(): momentum value must be a number')
            return
         end
         momentum = Value;
         if momentum > 1.0 || momentum < 0
          fprintf('runica(): momentum value is out of bounds [0,1]')
          return
         end
      elseif strcmp(Keyword,'sphering') || strcmp(Keyword,'sphereing') ...
                || strcmp(Keyword,'sphere')
         if ~ischar(Value)
           fprintf('runica(): sphering value must be on, off, or none')
           return
         else
           Value = lower(Value);
           if ~strcmp(Value,'on') && ~strcmp(Value,'off') && ~strcmp(Value,'none'),
             fprintf('runica(): sphering value must be on or off')
             return
           end
           sphering = Value;
         end
      elseif strcmp(Keyword,'bias')
         if ~ischar(Value)
           fprintf('runica(): bias value must be on or off')
           return
         else
           Value = lower(Value);
           if strcmp(Value,'on')
              biasflag = 1;
           elseif strcmp(Value,'off'),
              biasflag = 0;
           else
              fprintf('runica(): bias value must be on or off')
              return
           end
         end
      elseif strcmp(Keyword,'specgram') || strcmp(Keyword,'spec')

         if ~exist('specgram') < 2 % if ~exist or defined workspace variable
           fprintf(...
   'runica(): MATLAB Sig. Proc. Toolbox function "specgram" not found.\n')
           return
         end
         if ischar(Value)
           fprintf('runica(): specgram argument must be a vector')
           return
         end
         srate = Value(1);
         if (srate < 0)
             fprintf('runica(): specgram srate (%4.1f) must be >=0',srate)
             return
           end
         if length(Value)>1
           loHz = Value(2);
           if (loHz < 0 || loHz > srate/2)
             fprintf('runica(): specgram loHz must be >=0 and <= srate/2 (%4.1f)',srate/2)
             return
           end
         else
           loHz = 0; % default
         end
         if length(Value)>2
           hiHz = Value(3);
           if (hiHz < loHz || hiHz > srate/2)
             fprintf('runica(): specgram hiHz must be >=loHz (%4.1f) and <= srate/2 (%4.1f)',loHz,srate/2)
             return
           end
         else
           hiHz = srate/2; % default
         end
         if length(Value)>3
           Hzframes = Value(5);
           if (Hzframes<0 || Hzframes > size(data,2))
             fprintf('runica(): specgram frames must be >=0 and <= data length (%d)',size(data,2))
             return
           end
         else
           Hzframes = size(data,2); % default
         end
         if length(Value)>4
           Hzwinlen = Value(4);
           if rem(Hzframes,Hzwinlen) % if winlen doesn't divide frames
             fprintf('runica(): specgram Hzinc must divide frames (%d)',Hzframes)
             return
           end
         else
           Hzwinlen = Hzframes; % default
         end
         Specgramflag = 1; % set flag to perform specgram()

      elseif strcmp(Keyword,'extended') || strcmp(Keyword,'extend')
         if ischar(Value)
           fprintf('runica(): extended value must be an integer (+/-)')
           return
         else
           extended = 1;      % turn on extended-ICA
           extblocks = fix(Value); % number of blocks per kurt() compute
           if extblocks < 0
                nsub = -1*fix(extblocks);  % fix this many sub-Gauss comps
           elseif ~extblocks
                extended = 0;             % turn extended-ICA off
           elseif kurtsize>frames   % length of kurtosis calculation
                kurtsize = frames;
                if kurtsize < MIN_KURTSIZE
                   fprintf(...
   'runica() warning: kurtosis values inexact for << %d points.\n',...
                                                         MIN_KURTSIZE);
                end
           end
         end
      elseif strcmp(Keyword,'verbose')
         if ~ischar(Value)
            fprintf('runica(): verbose flag value must be on or off')
            return
         elseif strcmp(Value,'on')
             verbose = 1;
         elseif strcmp(Value,'off')
             verbose = 0;
         else
             fprintf('runica(): verbose flag value must be on or off')
             return
         end
      elseif strcmp(Keyword,'rndreset')
         if ischar(Value)
           if strcmp(Value,'yes')
             reset_randomseed = true;
           elseif strcmp(Value,'no')
             reset_randomseed = false;
           else
             fprintf('runica(): not using the reset_randomseed flag, it should be ''yes'',''no'',0, or 1');
           end
         else
           reset_randomseed = Value;
         end
      else
         fprintf('runica(): unknown flag')
         return
      end
   end

%
%%%%%%%%%%%%%%%%%%%%%%%% Initialize weights, etc. %%%%%%%%%%%%%%%%%%%%%%%%
%
if ~annealstep
  if ~extended
    annealstep = DEFAULT_ANNEALSTEP;     % defaults defined above
  else
    annealstep = DEFAULT_EXTANNEAL;       % defaults defined above
  end
end % else use annealstep from commandline

if ~annealdeg
    annealdeg  = DEFAULT_ANNEALDEG - momentum*90; % heuristic
    if annealdeg < 0
        annealdeg = 0;
    end
end
if ncomps >  chans || ncomps < 1
    fprintf('runica(): number of components must be 1 to %d.\n',chans);
    return
end
%
%%%%%%%%%%%%%%%%%%%%% Check keyword values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if frames<chans
    fprintf('runica(): data length (%d) < data channels (%d)!\n',frames,chans)
    return
elseif block < 2
    fprintf('runica(): block size %d too small!\n',block)
    return
elseif block > frames
    fprintf('runica(): block size exceeds data length!\n');
    return
elseif floor(epochs) ~= epochs
    fprintf('runica(): data length is not a multiple of the epoch length!\n');
    return
elseif nsub > ncomps
    fprintf('runica(): there can be at most %d sub-Gaussian components!\n',ncomps);
    return
end

if ~isempty(logfile)
    fid = fopen(logfile, 'w');
    if fid == -1, error('Cannot open logfile for writing'); end
else
    fid = [];
end
verb = verbose;

if weights ~= 0                    % initialize weights
  % starting weights are being passed to runica() from the commandline
    if  chans>ncomps && weights ~=0
        [r,c]=size(weights);
        if r~=ncomps || c~=chans
            fprintf('runica(): weight matrix must have %d rows, %d columns.\n', ...
                    chans,ncomps);
            return;
        end
    end
    icaprintf(verb,fid,'Using starting weight matrix named in argument list ...\n');
end;

%
% adjust nochange if necessary
%
if isnan(nochange)
    if ncomps > 32
        nochange = 1E-7;
        nochangeupdated = 1; % for fprinting purposes
    else
        nochangeupdated = 1; % for fprinting purposes
        nochange = DEFAULT_STOP;
    end
else
    nochangeupdated = 0;
end

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Process the data %%%%%%%%%%%%%%%%%%%%%%%%%%
%
icaprintf(verb,fid,'\nInput data size [%d,%d] = %d channels, %d frames/n', ...
          chans,frames,chans,frames);

if strcmp(pcaflag,'on')
    icaprintf(verb,fid,'After PCA dimension reduction,\n  finding ');
else
    icaprintf(verb,fid,'Finding ');
end
if ~extended
    icaprintf(verb,fid,'%d ICA components using logistic ICA.\n',ncomps);
else % if extended
icaprintf(verb,fid,'%d ICA components using extended ICA.\n',ncomps);
if extblocks > 0
    icaprintf(verb,fid,'Kurtosis will be calculated initially every %d blocks using %d data points.\n',...
              extblocks, kurtsize);
else
    icaprintf(verb,fid,'Kurtosis will not be calculated. Exactly %d sub-Gaussian components assumed.\n',nsub);
end
end
icaprintf(verb,fid,'Decomposing %d frames per ICA weight ((%d)^2 = %d weights, %d frames)\n',...
          floor(frames/ncomps.^2),ncomps.^2,frames);
icaprintf(verb,fid,'Initial learning rate will be %g, block size %d.\n',...
          lrate,block);
if momentum>0,
    icaprintf(verb,fid,'Momentum will be %g.\n',momentum);
end
icaprintf(verb,fid,'Learning rate will be multiplied by %g whenever angledelta >= %g deg.\n', ...
          annealstep,annealdeg);

if nochangeupdated
    icaprintf(verb,fid,'More than 32 channels: default stopping weight change 1E-7\n');
end
icaprintf(verb,fid,'Training will end when wchange < %g or after %d steps.\n', nochange,maxsteps);
if biasflag,
    icaprintf(verb,fid,'Online bias adjustment will be used.\n');
else
    icaprintf(verb,fid,'Online bias adjustment will not be used.\n');
end

%
%%%%%%%%%%%%%%%%% Remove overall row means of data %%%%%%%%%%%%%%%%%%%%%%%
%
icaprintf(verb,fid,'Removing mean of each channel ...\n');

%BLGBLGBLG replaced
% rowmeans = mean(data');
% data = data - rowmeans'*ones(1,frames);      % subtract row means
%BLGBLGBLG replacement starts
rowmeans = mean(data,2)'; %BLG
% data = data - rowmeans'*ones(1,frames);      % subtract row means
for iii=1:size(data,1) %avoids memory errors BLG
    data(iii,:)=data(iii,:)-rowmeans(iii);
end
%BLGBLGBLG replacement ends

icaprintf(verb,fid,'Final training data range: %g to %g\n', min(min(data)),max(max(data)));

%
%%%%%%%%%%%%%%%%%%% Perform PCA reduction %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
if strcmp(pcaflag,'on')
    icaprintf(verb,fid,'Reducing the data to %d principal dimensions...\n',ncomps);

    %BLGBLGBLG replaced
    %[eigenvectors,eigenvalues,data] = pcsquash(data,ncomps);
    % make data its projection onto the ncomps-dim principal subspace
    %BLGBLGBLG replacement starts
    %[eigenvectors,eigenvalues,data] = pcsquash(data,ncomps);
    % no need to re-subtract row-means, it was done a few lines above!
    PCdat2 = data';                    % transpose data
    [PCn,PCp]=size(PCdat2);                  % now p chans,n time points
    PCdat2=PCdat2/PCn;
    PCout=data*PCdat2;
    clear PCdat2;

    [PCV,PCD] = eig(PCout);                  % get eigenvectors/eigenvalues
    [PCeigenval,PCindex] = sort(diag(PCD));
    PCindex=rot90(rot90(PCindex));
    PCEigenValues=rot90(rot90(PCeigenval))';
    PCEigenVectors=PCV(:,PCindex);
    %PCCompressed = PCEigenVectors(:,1:ncomps)'*data;
    data = PCEigenVectors(:,1:ncomps)'*data;

    eigenvectors=PCEigenVectors;
    eigenvalues=PCEigenValues; %#ok<NASGU>

    clear PCn PCp PCout PCV PCD PCeigenval PCindex PCEigenValues PCEigenVectors
    %BLGBLGBLG replacement ends

end

%
%%%%%%%%%%%%%%%%%%% Perform specgram transformation %%%%%%%%%%%%%%%%%%%%%%%
%
if exist('Specgramflag') == 1
  % [P F T] = SPECGRAM(A,NFFT,Fs,WINDOW,NOVERLAP) % MATLAB Sig Proc Toolbox
  % Hzwinlen =  fix(srate/Hzinc); % CHANGED FROM THIS 12/18/00 -sm

  Hzfftlen = 2^(ceil(log(Hzwinlen)/log(2)));   % make FFT length next higher 2^k
  Hzoverlap = 0; % use sequential windows
  %
  % Get freqs and times from 1st channel analysis
  %
  [tmp,freqs,tms] = specgram(data(1,:),Hzfftlen,srate,Hzwinlen,Hzoverlap);

  fs = find(freqs>=loHz & freqs <= hiHz);
  icaprintf(verb,fid,'runica(): specified frequency range too narrow, exiting!\n');

  specdata = reshape(tmp(fs,:),1,length(fs)*size(tmp,2));
  specdata = [real(specdata) imag(specdata)];
     % fprintf('   size(fs) = %d,%d\n',size(fs,1),size(fs,2));
     % fprintf('   size(tmp) = %d,%d\n',size(tmp,1),size(tmp,2));
  %
  % Loop through remaining channels
  %
  for ch=2:chans
      [tmp] = specgram(data(ch,:),Hzwinlen,srate,Hzwinlen,Hzoverlap);
      tmp = reshape((tmp(fs,:)),1,length(fs)*size(tmp,2));
      specdata = [specdata;[real(tmp) imag(tmp)]]; % channels are rows
  end
  %
  % Print specgram confirmation and details
  %
  icaprintf(verb,fid,'Converted data to %d channels by %d=2*%dx%d points spectrogram data.\n',...
            chans,2*length(fs)*length(tms),length(fs),length(tms));
  if length(fs) > 1
      icaprintf(verb,fid,'   Low Hz %g, high Hz %g, Hz incr %g, window length %d\n',freqs(fs(1)),freqs(fs(end)),freqs(fs(2))-freqs(fs(1)),Hzwinlen);
  else
      icaprintf(verb,fid,'   Low Hz %g, high Hz %g, window length %d\n',freqs(fs(1)),freqs(fs(end)),Hzwinlen);
  end
  %
  % Replace data with specdata
  %
  data = specdata;
  datalength=size(data,2);
end
%
%%%%%%%%%%%%%%%%%%% Perform sphering %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if strcmp(sphering,'on'), %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  icaprintf(verb,fid,'Computing the sphering matrix...\n');
  sphere = 2.0*inv(sqrtm(double(cov(data')))); % find the "sphering" matrix = spher()
  if ~weights,
      icaprintf(verb,fid,'Starting weights are the identity matrix ...\n');
      weights = eye(ncomps,chans); % begin with the identity matrix
  else % weights given on commandline
      icaprintf(verb,fid,'Using starting weights named on commandline ...\n');
  end
  icaprintf(verb,fid,'Sphering the data ...\n');
  data = sphere*data; % decorrelate the electrode signals by 'sphereing' them

elseif strcmp(sphering,'off') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if ~weights % is starting weights not given
      icaprintf(verb,fid,'Using the sphering matrix as the starting weight matrix ...\n');
      icaprintf(verb,fid,'Returning the identity matrix in variable "sphere" ...\n');
      sphere = 2.0*inv(sqrtm(cov(data'))); % find the "sphering" matrix = spher()
      weights = eye(ncomps,chans)*sphere;  % begin with the identity matrix
      sphere = eye(chans);                 % return the identity matrix
  else % weights ~= 0
      icaprintf(verb,fid,'Using starting weights from commandline ...\n');
      icaprintf(verb,fid,'Returning the identity matrix in variable "sphere" ...\n');
      sphere = eye(chans);                 % return the identity matrix
  end
elseif strcmp(sphering,'none')
  sphere = eye(chans,chans);% return the identity matrix
  if ~weights
      icaprintf(verb,fid,'Starting weights are the identity matrix ...\n');
      icaprintf(verb,fid,'Returning the identity matrix in variable "sphere" ...\n');
      weights = eye(ncomps,chans); % begin with the identity matrix
  else % weights ~= 0
      icaprintf(verb,fid,'Using starting weights named on commandline ...\n');
      icaprintf(verb,fid,'Returning the identity matrix in variable "sphere" ...\n');
  end
  icaprintf(verb,fid,'Returned variable "sphere" will be the identity matrix.\n');
end
%
%%%%%%%%%%%%%%%%%%%%%%%% Initialize ICA training %%%%%%%%%%%%%%%%%%%%%%%%%
%
lastt=fix((datalength/block-1)*block+1);
BI=block*eye(ncomps,ncomps);
delta=zeros(1,chans*ncomps);
changes = [];
degconst = 180./pi;
startweights = weights;
prevweights = startweights;
oldweights = startweights;
prevwtchange = zeros(chans,ncomps);
oldwtchange = zeros(chans,ncomps);
lrates = zeros(1,maxsteps);
onesrow = ones(1,block);
bias = zeros(ncomps,1);
signs = ones(1,ncomps);    % initialize signs to nsub -1, rest +1
for k=1:nsub
    signs(k) = -1;
end
if extended && extblocks < 0,
    icaprintf(verb,fid,'Fixed extended-ICA sign assignments:  ');
    for k=1:ncomps
        icaprintf(verb,fid,'%d ',signs(k));
    end; icaprintf(verb,fid,'\n');
end
signs = diag(signs); % make a diagonal matrix
oldsigns = zeros(size(signs));
signcount = 0;              % counter for same-signs
signcounts = [];
urextblocks = extblocks;    % original value, for resets
old_kk = zeros(1,ncomps);   % for kurtosis momemtum

%
%%%%%%%% ICA training loop using the logistic sigmoid %%%%%%%%%%%%%%%%%%%
%
icaprintf(verb,fid,'Beginning ICA training ...');
if extended,
    icaprintf(verb,fid,' first training step may be slow ...\n');
else
    icaprintf(verb,fid,'\n');
end
step=0;
laststep=0;
blockno = 1;  % running block counter for kurtosis interrupts

warning('off', 'MATLAB:RandStream:ActivatingLegacyGenerators')
if reset_randomseed
    rand('state',sum(100*clock));  % set the random number generator state to
else
    rand('state', 0);
end                                % a position dependent on the system clock
warning('on', 'MATLAB:RandStream:ActivatingLegacyGenerators')

% interrupt figure
% ---------------
if strcmpi(interrupt, 'on')
    fig = figure('visible', 'off');
    supergui( fig, {1 1}, [], {'style' 'text' 'string' 'Press button to interrupt runica()' }, ...
              {'style' 'pushbutton' 'string' 'Interrupt' 'callback' 'setappdata(gcf, ''run'', 0);' } );
    set(fig, 'visible', 'on');
    setappdata(gcf, 'run', 1);

    if strcmpi(interrupt, 'on')
        drawnow;
    end
end

%% Compute ICA Weights
if biasflag && extended
    while step < maxsteps, %%% ICA step = pass through all the data %%%%%%%%%
        timeperm=randperm(datalength); % shuffle data order at each step

        for t=1:block:lastt, %%%%%%%%% ICA Training Block %%%%%%%%%%%%%%%%%%%
            if strcmpi(interrupt, 'on')
                drawnow;
                flag = getappdata(fig, 'run');
                if ~flag,
                    if ~isempty(fid), fclose(fid); end
                    close; error('USER ABORT');
                end
            end

            %% promote data block (only) to double to keep u and weights double
            u=weights*double(data(:,timeperm(t:t+block-1))) + bias*onesrow;

            y=tanh(u);
            weights = weights + lrate*(BI-signs*y*u'-u*u')*weights;
            bias = bias + lrate*sum((-2*y)')';  % for tanh() nonlin.

            if momentum > 0 %%%%%%%%% Add momentum %%%%%%%%%%%%%%%%%%%%%%%%%%%%
                weights = weights + momentum*prevwtchange;
                prevwtchange = weights-prevweights;
                prevweights = weights;
            end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            if max(max(abs(weights))) > MAX_WEIGHT
                wts_blowup = 1;
                change = nochange;
            end
            if ~wts_blowup
                %
                %%%%%%%%%%% Extended-ICA kurtosis estimation %%%%%%%%%%%%%%%%%%%%%
                %while step < maxsteps
                if extblocks > 0 && rem(blockno,extblocks) == 0,
                    % recompute signs vector using kurtosis
                    if kurtsize < frames % 12-22-99 rand() size suggestion by M. Spratling
                        rp = fix(rand(1,kurtsize)*datalength);  % pick random subset
                        % Accout for the possibility of a 0 generation by rand
                        ou = find(rp == 0);
                        while ~isempty(ou) % 1-11-00 suggestion by J. Foucher
                            rp(ou) = fix(rand(1,length(ou))*datalength);
                            ou = find(rp == 0);
                        end
                        partact=weights*double(data(:,rp(1:kurtsize)));
                    else                                        % for small data sets,
                        partact=weights*double(data);           % use whole data
                    end
                    m2=mean(partact'.^2).^2;
                    m4= mean(partact'.^4);
                    kk= (m4./m2)-3.0;                           % kurtosis estimates
                    if extmomentum
                        kk = extmomentum*old_kk + (1.0-extmomentum)*kk; % use momentum
                        old_kk = kk;
                    end
                    signs=diag(sign(kk+signsbias));             % pick component signs
                    if signs == oldsigns,
                        signcount = signcount+1;
                    else
                        signcount = 0;
                    end
                    oldsigns = signs;
                    signcounts = [signcounts signcount];
                    if signcount >= SIGNCOUNT_THRESHOLD,
                        extblocks = fix(extblocks * SIGNCOUNT_STEP);% make kurt() estimation
                        signcount = 0;                             % less frequent if sign
                    end                                         % is not changing
                end % extblocks > 0 & . . .
            end % if extended & ~wts_blowup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            blockno = blockno + 1;
            if wts_blowup
                break
            end
        end % for t=1:block:lastt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        if ~wts_blowup
            oldwtchange = weights-oldweights;
            step=step+1;
            %
            %%%%%%% Compute and print weight and update angle changes %%%%%%%%%
            %
            lrates(1,step) = lrate;
            angledelta=0.;
            delta=reshape(oldwtchange,1,chans*ncomps);
            change=delta*delta';
        end
        %
        %%%%%%%%%%%%%%%%%%%%%% Restart if weights blow up %%%%%%%%%%%%%%%%%%%%
        %
        if wts_blowup || isnan(change)|isinf(change),  % if weights blow up,
            icaprintf(verb,fid,'');
            step = 0;                          % start again
            change = nochange;
            wts_blowup = 0;                    % re-initialize variables
            blockno = 1;
            lrate = lrate*DEFAULT_RESTART_FAC; % with lower learning rate
            weights = startweights;            % and original weight matrix
            oldweights = startweights;
            change = nochange;
            oldwtchange = zeros(chans,ncomps);
            delta=zeros(1,chans*ncomps);
            olddelta = delta;
            extblocks = urextblocks;
            prevweights = startweights;
            prevwtchange = zeros(chans,ncomps);
            lrates = zeros(1,maxsteps);
            bias = zeros(ncomps,1);

            signs = ones(1,ncomps);    % initialize signs to nsub -1, rest +1
            for k=1:nsub
                signs(k) = -1;
            end
            signs = diag(signs); % make a diagonal matrix
            oldsigns = zeros(size(signs));;

            if lrate> MIN_LRATE
                r = rank(data); % determine if data rank is too low
                if r<ncomps
                    icaprintf(verb,fid,'Data has rank %d. Cannot compute %d components.\n',...
                        r,ncomps);
                    return
                else
                    icaprintf(verb,fid,...
                        'Lowering learning rate to %g and starting again.\n',lrate);
                end
            else
                icaprintf(verb,fid, ...
                    'runica(): QUITTING - weight matrix may not be invertible!\n');
                return;
            end
        else % if weights in bounds
            %
            %%%%%%%%%%%%% Print weight update information %%%%%%%%%%%%%%%%%%%%%%
            %
            if step> 2
                angledelta=acos((delta*olddelta')/sqrt(change*oldchange));
            end
            places = -floor(log10(nochange));
            icaprintf(verb,fid,'step %d - lrate %5f, wchange %8.8f, angledelta %4.1f deg\n', ...
                                step,      lrate,     change, degconst*angledelta);
            %
            %%%%%%%%%%%%%%%%%%%% Save current values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            changes = [changes change];
            oldweights = weights;
            %
            %%%%%%%%%%%%%%%%%%%% Anneal learning rate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            if degconst*angledelta > annealdeg,
                lrate = lrate*annealstep;          % anneal learning rate
                olddelta   = delta;                % accumulate angledelta until
                oldchange  = change;               %  annealdeg is reached
            elseif step == 1                     % on first step only
                olddelta   = delta;                % initialize
                oldchange  = change;
            end
            %
            %%%%%%%%%%%%%%%%%%%% Apply stopping rule %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            if step >2 && change < nochange,      % apply stopping rule
                laststep=step;
                step=maxsteps;                  % stop when weights stabilize
            elseif change > DEFAULT_BLOWUP,      % if weights blow up,
                lrate=lrate*DEFAULT_BLOWUP_FAC;    % keep trying
            end;                                 % with a smaller learning rate
        end; % end if weights in bounds

    end; % end while step < maxsteps (ICA Training) %%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% Compute ICA Weights
if biasflag && ~extended
    while step < maxsteps, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        timeperm=randperm(datalength); % shuffle data order at each step

        for t=1:block:lastt, %%%%%%%%% ICA Training Block %%%%%%%%%%%%%%%%%%%
            if strcmpi(interrupt, 'on')
                drawnow;
                flag = getappdata(fig, 'run');
                if ~flag,
                    if ~isempty(fid), fclose(fid); end
                    close; error('USER ABORT');
                end
            end

            u=weights*double(data(:,timeperm(t:t+block-1))) + bias*onesrow;
            y=1./(1+exp(-u));
            weights = weights + lrate*(BI+(1-2*y)*u')*weights;
            bias = bias + lrate*sum((1-2*y)')'; % for logistic nonlin. %

            if momentum > 0 %%%%%%%%% Add momentum %%%%%%%%%%%%%%%%%%%%%%%%%%%%
                weights = weights + momentum*prevwtchange;
                prevwtchange = weights-prevweights;
                prevweights = weights;
            end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            if max(max(abs(weights))) > MAX_WEIGHT
                wts_blowup = 1;
                change = nochange;
            end
            blockno = blockno + 1;
            if wts_blowup
                break
            end
        end % for t=1:block:lastt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        if ~wts_blowup
            oldwtchange = weights-oldweights;
            step=step+1;
            %
            %%%%%%% Compute and print weight and update angle changes %%%%%%%%%
            %
            lrates(1,step) = lrate;
            angledelta=0.;
            delta=reshape(oldwtchange,1,chans*ncomps);
            change=delta*delta';
        end
        %
        %%%%%%%%%%%%%%%%%%%%%% Restart if weights blow up %%%%%%%%%%%%%%%%%%%%
        %
        if wts_blowup || isnan(change)|isinf(change),  % if weights blow up,
            icaprintf(verb,fid,'');
            step = 0;                          % start again
            change = nochange;
            wts_blowup = 0;                    % re-initialize variables
            blockno = 1;
            lrate = lrate*DEFAULT_RESTART_FAC; % with lower learning rate
            weights = startweights;            % and original weight matrix
            oldweights = startweights;
            change = nochange;
            oldwtchange = zeros(chans,ncomps);
            delta=zeros(1,chans*ncomps);
            olddelta = delta;
            extblocks = urextblocks;
            prevweights = startweights;
            prevwtchange = zeros(chans,ncomps);
            lrates = zeros(1,maxsteps);
            bias = zeros(ncomps,1);
            if lrate> MIN_LRATE
                r = rank(data); % determine if data rank is too low
                if r<ncomps
                    icaprintf(verb,fid,'Data has rank %d. Cannot compute %d components.\n',r,ncomps);
                    return
                else
                    icaprintf(verb,fid,'Lowering learning rate to %g and starting again.\n',lrate);
                end
            else
                icaprintf(verb,fid,'runica(): QUITTING - weight matrix may not be invertible!\n');
                return;
            end
        else % if weights in bounds
            %
            %%%%%%%%%%%%% Print weight update information %%%%%%%%%%%%%%%%%%%%%%
            %
            if step> 2
                angledelta=acos((delta*olddelta')/sqrt(change*oldchange));
            end
            places = -floor(log10(nochange));
            icaprintf(verb,fid,'step %d - lrate %5f, wchange %8.8f, angledelta %4.1f deg\n', ...
                                step,      lrate,     change, degconst*angledelta);
            %
            %%%%%%%%%%%%%%%%%%%% Save current values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            changes = [changes change];
            oldweights = weights;
            %
            %%%%%%%%%%%%%%%%%%%% Anneal learning rate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            if degconst*angledelta > annealdeg,
                lrate = lrate*annealstep;          % anneal learning rate
                olddelta   = delta;                % accumulate angledelta until
                oldchange  = change;               %  annealdeg is reached
            elseif step == 1                     % on first step only
                olddelta   = delta;                % initialize
                oldchange  = change;
            end
            %
            %%%%%%%%%%%%%%%%%%%% Apply stopping rule %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            if step >2 && change < nochange,      % apply stopping rule
                laststep=step;
                step=maxsteps;                  % stop when weights stabilize
            elseif change > DEFAULT_BLOWUP,      % if weights blow up,
                lrate=lrate*DEFAULT_BLOWUP_FAC;    % keep trying
            end;                                 % with a smaller learning rate
        end; % end if weights in bounds

    end; % end while step < maxsteps (ICA Training) %%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% Compute ICA Weights
if ~biasflag && extended
    while step < maxsteps, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        timeperm=randperm(datalength); % shuffle data order at each step through data

        for t=1:block:lastt, %%%%%%%%% ICA Training Block %%%%%%%%%%%%%%%%%%%
            if strcmpi(interrupt, 'on')
                drawnow;
                flag = getappdata(fig, 'run');
                if ~flag,
                    if ~isempty(fid), fclose(fid); end
                    close; error('USER ABORT');
                end
            end

            u=weights*double(data(:,timeperm(t:t+block-1))); % promote block to dbl
            y=tanh(u);                                                       %
            weights = weights + lrate*(BI-signs*y*u'-u*u')*weights;

            if momentum > 0 %%%%%%%%% Add momentum %%%%%%%%%%%%%%%%%%%%%%%%%%%%
                weights = weights + momentum*prevwtchange;
                prevwtchange = weights-prevweights;
                prevweights = weights;
            end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            if max(max(abs(weights))) > MAX_WEIGHT
                wts_blowup = 1;
                change = nochange;
            end
            if ~wts_blowup
                %
                %%%%%%%%%%% Extended-ICA kurtosis estimation %%%%%%%%%%%%%%%%%%%%%
                %while step < maxsteps
                if extblocks > 0 && rem(blockno,extblocks) == 0,
                    % recompute signs vector using kurtosis
                    if kurtsize < frames % 12-22-99 rand() size suggestion by M. Spratling
                        rp = fix(rand(1,kurtsize)*datalength);  % pick random subset
                        % Accout for the possibility of a 0 generation by rand
                        ou = find(rp == 0);
                        while ~isempty(ou) % 1-11-00 suggestion by J. Foucher
                            rp(ou) = fix(rand(1,length(ou))*datalength);
                            ou = find(rp == 0);
                        end
                        partact=weights*double(data(:,rp(1:kurtsize)));
                    else                                        % for small data sets,
                        partact=weights*double(data);           % use whole data
                    end
                    m2=mean(partact'.^2).^2;
                    m4= mean(partact'.^4);
                    kk= (m4./m2)-3.0;                           % kurtosis estimates
                    if extmomentum
                        kk = extmomentum*old_kk + (1.0-extmomentum)*kk; % use momentum
                        old_kk = kk;
                    end
                    signs=diag(sign(kk+signsbias));             % pick component signs
                    if signs == oldsigns,
                        signcount = signcount+1;
                    else
                        signcount = 0;
                    end
                    oldsigns = signs;
                    signcounts = [signcounts signcount];
                    if signcount >= SIGNCOUNT_THRESHOLD,
                        extblocks = fix(extblocks * SIGNCOUNT_STEP);% make kurt() estimation
                        signcount = 0;                             % less frequent if sign
                    end                                         % is not changing
                end % extblocks > 0 & . . .
            end % if ~wts_blowup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            blockno = blockno + 1;
            if wts_blowup
                break
            end
        end % for t=1:block:lastt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        if ~wts_blowup
            oldwtchange = weights-oldweights;
            step=step+1;
            %
            %%%%%%% Compute and print weight and update angle changes %%%%%%%%%
            %
            lrates(1,step) = lrate;
            angledelta=0.;
            delta=reshape(oldwtchange,1,chans*ncomps);
            change=delta*delta';
        end
        %
        %%%%%%%%%%%%%%%%%%%%%% Restart if weights blow up %%%%%%%%%%%%%%%%%%%%
        %
        if wts_blowup || isnan(change)|isinf(change),  % if weights blow up,
            icaprintf(verb,fid,'');
            step = 0;                          % start again
            change = nochange;
            wts_blowup = 0;                    % re-initialize variables
            blockno = 1;
            lrate = lrate*DEFAULT_RESTART_FAC; % with lower learning rate
            weights = startweights;            % and original weight matrix
            oldweights = startweights;
            change = nochange;
            oldwtchange = zeros(chans,ncomps);
            delta=zeros(1,chans*ncomps);
            olddelta = delta;
            extblocks = urextblocks;
            prevweights = startweights;
            prevwtchange = zeros(chans,ncomps);
            lrates = zeros(1,maxsteps);
            bias = zeros(ncomps,1);
            signs = ones(1,ncomps);    % initialize signs to nsub -1, rest +1
            for k=1:nsub
                signs(k) = -1;
            end
            signs = diag(signs); % make a diagonal matrix
            oldsigns = zeros(size(signs));
            if lrate> MIN_LRATE
                r = rank(data); % find whether data rank is too low
                if r<ncomps
                    icaprintf(verb,fid,'Data has rank %d. Cannot compute %d components.\n',...
                        r,ncomps);
                    return
                else
                    icaprintf(verb,fid,...
                        'Lowering learning rate to %g and starting again.\n',lrate);
                end
            else
                icaprintf(verb,fid, ...
                    'runica(): QUITTING - weight matrix may not be invertible!\n');
                return;
            end
        else % if weights in bounds
            %
            %%%%%%%%%%%%% Print weight update information %%%%%%%%%%%%%%%%%%%%%%
            %
            if step> 2
                angledelta=acos((delta*olddelta')/sqrt(change*oldchange));
            end
            places = -floor(log10(nochange));
            icaprintf(verb,fid,'step %d - lrate %5f, wchange %8.8f, angledelta %4.1f deg\n', ...
                                step,      lrate,     change, degconst*angledelta);
            %
            %%%%%%%%%%%%%%%%%%%% Save current values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            changes = [changes change];
            oldweights = weights;
            %
            %%%%%%%%%%%%%%%%%%%% Anneal learning rate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            if degconst*angledelta > annealdeg,
                lrate = lrate*annealstep;          % anneal learning rate
                olddelta   = delta;                % accumulate angledelta until
                oldchange  = change;               %  annealdeg is reached
            elseif step == 1                     % on first step only
                olddelta   = delta;                % initialize
                oldchange  = change;
            end
            %
            %%%%%%%%%%%%%%%%%%%% Apply stopping rule %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            if step >2 && change < nochange,      % apply stopping rule
                laststep=step;
                step=maxsteps;                  % stop when weights stabilize
            elseif change > DEFAULT_BLOWUP,      % if weights blow up,
                lrate=lrate*DEFAULT_BLOWUP_FAC;    % keep trying
            end;                                 % with a smaller learning rate
        end; % end if weights in bounds

    end; % end while step < maxsteps (ICA Training) %%%%%%%%%%%%%%%%%%%%%%%%%
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Compute ICA Weights
if ~biasflag && ~extended
    while step < maxsteps, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        timeperm=randperm(datalength); % shuffle data order at each step

        for t=1:block:lastt, %%%%%%%%% ICA Training Block %%%%%%%%%%%%%%%%%%%
            if strcmpi(interrupt, 'on')
                drawnow;
                flag = getappdata(fig, 'run');
                if ~flag,
                    if ~isempty(fid), fclose(fid); end
                    close; error('USER ABORT');
                end
            end
            u=weights*double(data(:,timeperm(t:t+block-1)));
            y=1./(1+exp(-u));                                                %
            weights = weights + lrate*(BI+(1-2*y)*u')*weights;

            if momentum > 0 %%%%%%%%% Add momentum %%%%%%%%%%%%%%%%%%%%%%%%%%%%
                weights = weights + momentum*prevwtchange;
                prevwtchange = weights-prevweights;
                prevweights = weights;
            end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            if max(max(abs(weights))) > MAX_WEIGHT
                wts_blowup = 1;
                change = nochange;
            end

            blockno = blockno + 1;
            if wts_blowup
                break
            end
        end % for t=1:block:lastt %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        if ~wts_blowup
            oldwtchange = weights-oldweights;
            step=step+1;
            %
            %%%%%%% Compute and print weight and update angle changes %%%%%%%%%
            %
            lrates(1,step) = lrate;
            angledelta=0.;
            delta=reshape(oldwtchange,1,chans*ncomps);
            change=delta*delta';
        end
        %
        %%%%%%%%%%%%%%%%%%%%%% Restart if weights blow up %%%%%%%%%%%%%%%%%%%%
        %
        if wts_blowup || isnan(change)|isinf(change),  % if weights blow up,
            icaprintf(verb,fid,'');
            step = 0;                          % start again
            change = nochange;
            wts_blowup = 0;                    % re-initialize variables
            blockno = 1;
            lrate = lrate*DEFAULT_RESTART_FAC; % with lower learning rate
            weights = startweights;            % and original weight matrix
            oldweights = startweights;
            change = nochange;
            oldwtchange = zeros(chans,ncomps);
            delta=zeros(1,chans*ncomps);
            olddelta = delta;
            extblocks = urextblocks;
            prevweights = startweights;
            prevwtchange = zeros(chans,ncomps);
            lrates = zeros(1,maxsteps);
            bias = zeros(ncomps,1);

            if lrate> MIN_LRATE
                r = rank(data); % find whether data rank is too low
                if r<ncomps
                    icaprintf(verb,fid,'Data has rank %d. Cannot compute %d components.\n',...
                        r,ncomps);
                    return
                else
                    icaprintf(verb,fid,...
                        'Lowering learning rate to %g and starting again.\n',lrate);
                end
            else
                icaprintf(verb,fid, ...
                    'runica(): QUITTING - weight matrix may not be invertible!\n');
                return;
            end
        else % if weights in bounds
            %
            %%%%%%%%%%%%% Print weight update information %%%%%%%%%%%%%%%%%%%%%%
            %
            if step> 2
                angledelta=acos((delta*olddelta')/sqrt(change*oldchange));
            end
            places = -floor(log10(nochange));
            icaprintf(verb,fid,'step %d - lrate %5f, wchange %8.8f, angledelta %4.1f deg\n', ...
                                step,      lrate,     change, degconst*angledelta);
            %
            %%%%%%%%%%%%%%%%%%%% Save current values %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            changes = [changes change];
            oldweights = weights;
            %
            %%%%%%%%%%%%%%%%%%%% Anneal learning rate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            if degconst*angledelta > annealdeg,
                lrate = lrate*annealstep;          % anneal learning rate
                olddelta   = delta;                % accumulate angledelta until
                oldchange  = change;               %  annealdeg is reached
            elseif step == 1                     % on first step only
                olddelta   = delta;                % initialize
                oldchange  = change;
            end
            %
            %%%%%%%%%%%%%%%%%%%% Apply stopping rule %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            if step >2 && change < nochange,      % apply stopping rule
                laststep=step;
                step=maxsteps;                  % stop when weights stabilize
            elseif change > DEFAULT_BLOWUP,      % if weights blow up,
                lrate=lrate*DEFAULT_BLOWUP_FAC;    % keep trying
            end;                                 % with a smaller learning rate
        end; % end if weights in bounds

    end; % end while step < maxsteps (ICA Training) %%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% Finalize Computed Data for Output

if strcmpi(interrupt, 'on')
    close(fig);
end

  if ~laststep
    laststep = step;
  end
  lrates = lrates(1,1:laststep);           % truncate lrate history vector

  %
  %%%%%%%%%%%%%% Orient components towards max positive activation %%%%%%
  %
  if nargout > 6 || strcmp(posactflag,'on')
      % make activations from sphered and pca'd data; -sm 7/05
      % add back the row means removed from data before sphering
      if strcmp(pcaflag,'off')
          sr = sphere * rowmeans';
          for r = 1:ncomps
              data(r,:) = data(r,:)+sr(r); % add back row means
          end
          data = weights*data; % OK in single
      else
          ser = sphere*eigenvectors(:,1:ncomps)'*rowmeans';
          for r = 1:ncomps
              data(r,:) = data(r,:)+ser(r); % add back row means
          end
          data = weights*data; % OK in single
      end
  end
  %
  % NOTE: Now 'data' are the component activations = weights*sphere*raw_data
  %
  %
  %%%%%%%%%%%%%% If pcaflag, compose PCA and ICA matrices %%%%%%%%%%%%%%%
  %
  if strcmp(pcaflag,'on')
    icaprintf(verb,fid,'Composing the eigenvector, weights, and sphere matrices\n');
    icaprintf(verb,fid,'  into a single rectangular weights matrix; sphere=eye(%d)\n'...
                                                                  ,chans);
    weights= weights*sphere*eigenvectors(:,1:ncomps)';
    sphere = eye(urchans);
  end
  %
  %%%%%% Sort components in descending order of max projected variance %%%%
  %
  icaprintf(verb,fid,'Sorting components in descending order of mean projected variance ...\n');
  %
  %%%%%%%%%%%%%%%%%%%% Find mean variances %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % meanvar  = zeros(ncomps,1);      % size of the projections
  if ncomps == urchans % if weights are square . . .
      winv = inv(weights*sphere);
  else
      icaprintf(verb,fid,'Using pseudo-inverse of weight matrix to rank order component projections.\n');
      winv = pinv(weights*sphere);
  end
  %
  % compute variances without backprojecting to save time and memory -sm 7/05
  %
  meanvar = sum(winv.^2).*sum((data').^2)/((chans*frames)-1); % from Rey Ramirez 8/07
  %
  %%%%%%%%%%%%%% Sort components by mean variance %%%%%%%%%%%%%%%%%%%%%%%%
  %
  [sortvar, windex] = sort(meanvar);
  windex = windex(ncomps:-1:1); % order large to small
  meanvar = meanvar(windex);
  %
  %%%%%%%%%%%% re-orient max(abs(activations)) to >=0 ('posact') %%%%%%%%
  %
  if strcmp(posactflag,'on') % default is now off to save processing and memory
      icaprintf(verb,fid,'Making the max(abs(activations)) positive ...\n');
      [tmp ix] = max(abs(data')); % = max abs activations
      signsflipped = 0;
      for r=1:ncomps
         if sign(data(r,ix(r))) < 0
            if nargout>6  % if activations are to be returned (only)
               data(r,:) = -1*data(r,:);  % flip activations so max(abs()) is >= 0
            end
            winv(:,r) = -1*winv(:,r);  % flip component maps
            signsflipped = 1;
         end
      end
      if signsflipped == 1
          weights = pinv(winv)*inv(sphere); % re-invert the component maps
      end

      % [data,winvout,weights] = posact(data,weights); % overwrite data with activations
      % changes signs of activations (now = data) and weights
      % to make activations (data) net rms-positive
      % can call this outside of runica() - though it is inefficient!
  end
  %
  %%%%%%%%%%%%%%%%%%%%% Filter data using final weights %%%%%%%%%%%%%%%%%%
  %
  if nargout>6, % if activations are to be returned
      icaprintf(verb,fid,'Permuting the activation wave forms ...\n');
      data = data(windex,:); % data is now activations -sm 7/05
  else
      clear data
  end
  weights = weights(windex,:);% reorder the weight matrix
  bias  = bias(windex);       % reorder them
  signs = diag(signs);        % vectorize the signs matrix
  signs = signs(windex);      % reorder them

  if ~isempty(fid), fclose(fid); end; % close logfile

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

return

% printing functions
% ------------------
function icaprintf(verb,fid, varargin);
    if verb
        if ~isempty(fid)
            fprintf(fid, varargin{:});
        end;
        fprintf(varargin{:});
    end


Officers responded to a family fight involving three women and two men at a residence on Crestview Drive at approximately 8 p.m. Nov. 23, according to a probable cause statement filed by the Santa Clara-Ivins Police Department in support of the arrest. Six children were also present in the home at the time.

The argument and subsequent physical altercation triggered a female pitbull to begin attacking and biting individuals in the home, including a 38-year-old man, a 31-year-old woman and two children, ages 13 and 3, the report states.

The pitbull had been caring for a new litter of puppies in a metal dog crate inside the home when she became agitated, broke out of the crate and began the attack, resulting in the four individuals being transported by ambulance to Dixie Regional Medical Center in St. George to be treated for their injuries.

Two additional family members rode along in the ambulances but were not injured, Santa Clara-Ivins Police Public Information Officer Randy Hancey said Friday, clarifying earlier statements that six injured were transported.

Animal control officers removed an adult male pitbull, the female pitbull and her litter of pups from the home, Hancey said, adding that he did not know if the male pitbull was involved in the attack. He said the pitbulls would be quarantined for a period of time as a matter of standard protocol.

Management personnel at the Ivins Animal Shelter declined both comment and access to the dogs, referring St. George News to the Police Department.

A 31-year-old woman involved in the fight has not been named in this report as she has not been arrested or charged as this report publishes. The 31-year-old woman and Navarro had previously lived at the Crestview Drive residence with Lopez, Montez and Campos, who currently live in the home.

They had allegedly all been consuming large amounts of alcohol before the 31-year-old woman got into a fight with Montez and Campos, the arresting officer wrote in the sworn statement, but the statement did not specify what the dispute was about.

When officers first arrived on scene, they saw Montez and the 31-year-old woman in the living room of the residence punching each other with closed fists, the report states.

At some point, prior to officers arriving on scene, Lopez got involved in the fight, punched Navarro in the face and threatened Navarro’s life by stating he “was going to kill him,” according to the statement.

“I also found Isaac to admit to not remembering much of the incident due to heavy alcohol consumption,” the officer wrote, adding that Montez and Lopez also said all six children were present during the altercation.

“While at the residence, officers detected the smell of marijuana,” the report states. “During questioning, (Lopez) admitted the marijuana to be his and gave officers custody of the marijuana in a silver thermos.”

On Wednesday, once officers had a chance to conduct interviews with everyone involved, Lopez, Navarro, Campos and Montez were arrested and booked into the Washington County Purgatory Correctional Facility.

The Washington County Attorney’s Office filed an Information in 5th District Court accusing Lopez of three class A misdemeanor counts of drug possession, threat of violence and assault, along with six additional class A misdemeanor counts of domestic violence in the presence of a child, class B misdemeanor possession of drug paraphernalia and class C misdemeanor intoxication.

Lopez has since been released from custody on $16,790 bail pending trial. His arraignment is scheduled for Dec. 1 before Judge G. Michael Westfall.

According to court documents, since 2013, Lopez has been arrested in Utah on a variety of charges including disorderly conduct, criminal mischief for intentionally destroying property, drug-related charges and domestic violence in the presence of a child.

The county attorney’s office has charged Navarro with class A misdemeanor assault, six class A misdemeanor counts of domestic violence in the presence of a child, and class C misdemeanor intoxication.

Navarro was also found to have three active warrants for his arrest issued by the Washington County 5th District Court relating to driving without a valid driver’s license.

Campos has been charged with class B misdemeanor assault and six class B misdemeanor counts of domestic violence in the presence of a child.

Campos appeared in Santa Clara Justice Court before Judge Ken Armstrong Thursday afternoon and pleaded not guilty to the charges. A pretrial conference for her case is scheduled for Jan. 26, 2018.

While booking information shows Montez was arrested on eight class B misdemeanor counts for interfering with arrest, negligent child abuse and domestic violence in the presence of a child, court records do not indicate charges had officially been filed as this report is published.

In a preliminary interview with Linda Montez she told St. George News she will defend the charges against her as her involvement was breaking up a fight between other people in her own home.

Montez also said that the two dogs, Shadow and Nala, have no previous incidents of aggression.

This is a developing story.

This report is based on preliminary information provided by law enforcement or other emergency responders and may not contain the full scope of findings.

Free News Delivery by Email

Would you like to have the day's news stories delivered right to your inbox every evening? Enter your email below to start!

About the Author

Kimberly Scott is a lifetime resident of Southern Utah. In 2013, Kimberly joined St. George News as a vital member of its editing, reporting and administrative team. She is passionate about engaging communities through writing and is dedicated to providing complete and accurate coverage of both anticipated and breaking news.

Just for anyone info we are not from south of the border ok I’m a US citizen just like all you folks that like to run your mouths! So the story is lies that’s is not what happens at all ok I can say what happens but I’ll let my lawyer talk so for any ones comments Linda Montez, Isaac Lopez and Griselda Campos we are innocent and we all have lawyers and are gonna prove our innocence so stfu and mind your buisness


Q:

Best way to duplicate Linux server

My company has an old CentOS 5 system that we need to duplicate for testing. I'm running dd on 2 of the partitions which are small and saving them to /home, but the /home partition can't save to itself or another partition as there isn't enough room.
/home is only about 7% full so copying all of its files to a tarball should be easy enough, but I'm wondering if there's a better way to do this.
I can't bring the machine offline and I can't hook up a USB drive externally and copy to files that way so I have to either use SCP, dd or tar most likely.
Any suggestions? Here's a quick breakdown of df -h
Filesystem  Size   Used  Left usage% Mount point
/dev/sda2    97G   17G   75G  19%     /
/dev/sda1    99M   15M   80M  16%     /boot
tmpfs       2.0G     0  2.0G   0%     /dev/shm
/dev/sda3   352G   22G  312G   7%     /home

A:

There are definitely "better" ways to do this, but how you go about doing this comes down to what resources you have available. I'll list a few different ways and you can decide which works best for you.
To begin, copying data in /home to a tarball is pointless. You can't save it locally, which means that to do this, it will require a remote server that you can send the data to. Preferably, a server on the same network as this computer that is accessible via SSH. Frankly, making tarballs of the other partitions and putting them in /home is a waste of time.
If you're looking to duplicate this server, what I would do is use dd to make an image of /dev/sda and send the data over SSH to another server. The two best ways of doing this are to either image another hard drive at the same time or save the image as a file.
dd if=/dev/sda bs=4096 conv=notrunc,noerror,sync | ssh user@machine "dd of=/dev/sdX"
This will clone the entire drive, including MBR (and therefore bootloader), all partitions, UUIDs, and data, sending it over SSH to another machine and imaging a drive of your choosing. The options I suggested do the following:

notrunc or 'do not truncate' maintains data integrity by instructing
dd not to truncate any data.
noerror instructs dd to continue operation, ignoring all read errors.
Default behavior for dd is to halt at any error.
sync writes zeroes for read errors, so data offsets stay in sync.
bs=4096 sets the block size to 4k, an optimal size for hard disk
read/write efficiency and therefore, cloning speed.

For this to work, the user you use log in with SSH as needs to either be root, or part of the disk group. /dev/sdaX should be mode 660, giving group read/write access. The default group owner should be disk.
Alternatively, you can create an image that you can use to write directly to the disk.
dd if=/dev/sda bs=4096 conv=notrunc,noerror,sync | ssh user@machine "dd of=/path/to/wherever/server.img"
To write the image to a hard drive from the computer it is saved on, you can use the following command.
dd if=/path/to/server.img bs=4096 conv=notrunc,noerror,sync of=/dev/sdX
Lastly, if space is a concern, you can gzip the image before sending it over SSH.
dd if=/dev/hda conv=sync,noerror bs=4096 | gzip -c | ssh user@machine "of=/path/to/wherever/server.img.gz"
and writing it by
gunzip -c /path/to/server.img.gz | dd of=/dev/sdX
Lastly, if you want to pull individual files from the image without writing it to another drive, you can mount it. However, this won't work if it's compressed.
You'll have to fdisk -lu /path/to/server.img to get the offset of the individual partitions. Once you do that, you need the starting sector of the partition you want to mount. For the sake of argument, let us say that the partition we want to access starts on sector 122441728.
$ fdisk -lu server.img

Disk server.img: 320.1 GB, 320072933376 bytes
255 heads, 63 sectors/track, 38913 cylinders, total 625142448 sectors
Units = sectors of 1 * 512 = 512 bytes
Sector size (logical/physical): 512 bytes / 512 bytes
I/O size (minimum/optimal): 512 bytes / 512 bytes
Disk identifier: 0x86308630

   Device Boot      Start         End      Blocks   Id  System
/dev/sda1   *          63   114441469    57220703+  83  Linux
/dev/sda2       114442240   122441727     3999744   82  Linux swap / Solaris
/dev/sda3       122441728   625141759   251350016   83  Linux

You would then use mount -o ro,loop,offset=$((122441728 * 512)) server.img /mnt to mount the third partition to /mnt.
Sources:

https://wiki.archlinux.org/index.php/Disk_Cloning#Create_disk_image
http://madduck.net/blog/2006.10.20:loop-mounting-partitions-from-a-disk-image/


[Elimination of succinate and acetate synthesis in recombinant Escherichia coli for D-lactate production].
When Escherichia coli CICIM B0013-030 (B0013, ack-pta, pps, pflB) was used for D-lactate production, succinate and acetate were the main byproducts (as much as 11.9 and 7.1% the amount of lactate respectively). In order to decrease the byproduct levels, we inactivated succinate and acetate synthesis in B0013-030. Two recombinant plasmids containing mutation cassettes of frdA::difGm and tdcDE::difGm respectively were constructed first. The mutation cassettes were used to delete the target genes on the chromosomal by Red recombination. Subsequently, the antibiotic resistance gene was excised from the chromosomal by Xer recombination. Thereby, mutants B0013-040B (B0013-030, frdA) and B0013-050B (B0013-040B, tdcDE) were produced. D-lactate producing abilities of the engineered strains were tested both in shake flasks and in bioreactors using two-phase fermentation (aerobic growth and anaerobic fermentation) with glucose as the sole carbon source. When fermentation was carried out in shake flasks, inactivation of frdA in B0013-030 to produce B0013-040B reduced succinate accumulation by 80.8%. When tested in a 7-liter bioreactor, B0013-040B accumulated 114.5 g/L D-lactate of over 99.9% optical purity. However, 1.0 g/L succinate and 5.4 g/L acetate still remained in the broth. Further inactivation of tdcD and tdcE genes in B0013-040B to produce B0013-050B decreased acetate and succinate accumulation to 0.4 g/L and 0.4 g/L respectively, and lactate titer was as much as 111.9 g/L (tested in the 7-liter bioreactor). In lightof the lower byproduct levels and high lactate production, strain B00 13-050B may prove useful for D-lactate production.


Here’s a case study in why the awakened man should never put himself into the compromised position of being married these days. This case is still winding its way through the vagina-worshiping court system, but it sounds a number of alarm bells to the man who would glean knowledge from its details.

The press is having a field day with the Emmanuel Sanders divorce. We will look at it from a Red Pill perspective, yet another cautionary tale of what putting your scrotum into a legal vice grip yields when your bitch gets ready to flush you down the toilet.

Here goes. (This is all just one man’s opinion, of course.)

Their second baby hadn’t even been born yet, and already Denver Broncos Wide Receiver Emmanuel Sanders’ wife Gabriella was positioning herself for the kill. She was ready to cash in on his ass as her divorce lawyers had already begun writing up boilerplate legal screeds for the gynocentric, Anglo-American kangaroo courts. Gabriella served Sanders with divorce papers even before she squeezed out her annuity… er… child in the delivery room.

In case you haven’t noticed, children are nothing more than meal tickets to women in modern day Anglo-America. Fathers and fathers’ rights be damned. If women think they can’t cash in on their children, they simply go kill them at the abortion clinic. But, being married to an NFL player obviously ensures a nice, consistent stream of child support cash enforced by legal indentured servitude. (Pay up or go to jail.) We are living smack dab in the middle of one of the most insane socioeconomic systems in the history of mankind.

How much does she want? As reported by Meninist:

Gabriella is looking to score big as she is asking the courts for custody of their two children, a restraining order against Sanders, and she wants a “disproportionate” share of their assets. Ironically, Sanders is scheduled to make $13.2 million (includes $10 Million dollar signing bonus) this year from his 3 year $33 million dollar contract.

Disproportionate share, eh? Tell me again how women are not prostitutes and marriage is not legalized prostitution.

Witness, once again White Knights and fembots, The Predatory Female in action. Gabriella timed the divorce perfectly. She filed for divorce right before she pushed out her 18-year annuity that Sanders fathered. But she doesn’t just want a fat child support check. Oh, no. She wants a fat slice of the bonus he no doubt busts his ass to earn.

Sanders even reportedly said he was “blindsided” by the divorce. Ha. Yep, women will do that. They always hide their true feelings and emotions from men until they go for the jugular. This should come as absolutely no surprise to a regular reader of TNMM. Cynical as the truth may be, the truth is the truth. Women are financial predators of men, and in a society that has gone completely insane and let them frivorce rape with impunity, they’re cleaning out any men stupid enough to marry them.

Laughably, she is citing his side chicks as a reason she should be able to cash in. (I want money for the sex I gave you rule applies again!) From celebrity gossip rag TMZ:

The estranged wife of Denver Broncos WR Emmanuel Sanders has declared WAR on the NFL star — claiming he lied, cheated and spent THOUSANDS on chicks he was banging on the side.

Umm…THOUSANDS to an NFL player is like dropping a $50 to the average working slob. It seems Sanders was another fool who soon parted with his money to pay for these girls rather than a PUA. A player on the field, and another Beta with more money than common sense when it comes to women. TMZ rages on:

It’s all in ANGRY court docs filed by Gabriella Sanders — who blasts the father of her 2 children with a metaphorical bazooka right from the get go. “When he’s not playing football on the field, he ‘plays the field’ with numerous women whom he is or has committed adultery.”

Big shocker. Knowing an NFL player will more often than not have multiple harems because women will lay down for rich and powerful men on command, nay throw themselves in front of his cock, she sets herself up as…you guessed it…the victim while she lends him enough rope to hang himself.

Gabby lays out multiple marital “atrocities” allegedly committed by her husband – and claims he spent “thousands upon thousands of dollars on girlfriends and wasting the community estate, even purchasing a vehicle for one of his illicit affairs.”

There is no doubt a personal investigator was used to follow Sanders around in order to gather information on him while she planned her frivorce rape lotto victory.

Incidents like these remind me of the Incident at the Bank I have written about numerous times, in which two Anglobitches compared notes about how much money their vaginas had rewarded them in divorce court as they cashed out at a Latin American teller. My opinion of women was already pretty low at that point. But I was forced to downgrade my opinion after witnessing that cackling display.

The long and short of it is only a fool marries an Anglobitch or a woman who has adapted to Anglo culture in today’s legal and social climate. That humble advice could have saved Sanders millions. Any other rich guys who I talk out of getting married, feel free to send me a nice contribution. I’ll take a fraction of what these hos will want out of you for my consulting fee.

The story of the Scorpion and the Frog applies: Men should not be asking why they were frivorce raped by a creature whose entire nature is built on dissimulation and acting. What they should be asking is why they were stupid enough to fall for the ruse. Taking advantage of men for material and economic gain is in women’s nature. And America was insane to let these instincts return, unchecked. It’s all just a downward spiral from here.

We’ll be on the sidelines keeping you advised of some of the more egregious frivorce rape lotto cases. These happen every day, it’s not like they’re uncommon.

Ergo, women always want money for sex. All the pleasantries around the exchange and myths of “love” and “marriage” and “commitment” are just sophistry and dissimulation. Find a divorce case in which women don’t want monetary remuneration for services performed, even in the YouGoGrrl age, and you will have found the proverbial needle in the haystack.

Help us grow by making a purchase from our Recommended Reading and Viewing page or our Politically Incorrect Apparel and Merchandise page or buy anything from Amazon using this link. You can also Sponsor The New Modern Man for as little as $1 a month. This article is the opinion of its author.


Editorial: The world is watching for quality

In its long battle for international recognition, Taiwan has no greater asset than its hard-won democracy. More than anything, its democratic achievement is what distinguishes it from the authoritarian system in China. It is therefore crucial that this democracy be seen by the international community to be a functioning one worth protecting.

With crackdowns on human-rights advocates, curtailed freedom of the press and threats of illegal use of force against another nation, Beijing continues to provide ample contrasts with Taiwan, and those differences are part of the means by which the latter can express its identity.

But the determination of one's identity through differentiation can only accomplish so much. A nation cannot define who and what it is solely by focusing on what it is not. In fact, to find its true self, it must also make a statement about what it is.

And what it is is much more than what differentiates it from China -- a democratic system versus an authoritarian regime -- and comes instead from how it has developed the gift of democracy.

As such, next year's elections must be about more than the fact that elections can be held in the first place. Beyond process alone, it is the substance and quality of the elections that constitute the true health check of the nation and that provide it with the opportunity to underscore its value to the world.

Hence the need for voters next year to be offered choices among contending politicians who have been cleared, through impartial investigations, of wrongdoing or corruption. Nothing could be more harmful to Taiwan's survival and ambitions for sovereignty than for the world to see that its people value their democracy so lightly as to bring a crook to power. Should the world come to the conclusion that Taiwanese opted to elect an individual who has openly lied about his powers while in office, it will no longer feel the obligation to protect that society from anti-democratic encroachment.

After all, we cannot expect the rest of the world to care about Taiwan's democracy if its very people are incapable of seeing to it that it remains healthy and truly representative.

Just as Taiwan's military allies need to see that it is willing to do what it must and acquire what it needs to defend its territory from military aggression, its diplomatic allies -- and in fact even those who side with Beijing -- need to see that Taiwanese are ready to do what they must to protect their precious democracy. Commitment to that principle -- the process itself as well as the quality of the choices given a people -- is just as important as commitment to national defense.

Canada recently showed what it means to protect the quality of a democratic system by pushing for the investigation of high-level government corruption. The process ultimately led to the demise of prime minister Paul Martin's Liberal government. Regardless of whether the new Conservative government was the best thing for Canada, the country was mature enough to rid itself of elements whose presence was harmful to its democracy.

Taiwan should follow that example and ensure that allegations of corruption involving the upper echelons of government are fully explored. The prosecutors' decision to appeal the ruling on former Taipei mayor Ma Ying-jeou (馬英九) in the special funds scandal is a step in that direction.


Wes Jarvis

Wesley Herbert Jarvis (born May 30, 1958 in Toronto, Ontario) is a former Canadian ice hockey player, and is a cousin of former NHLer Doug Jarvis.

Hockey career
As a youth, Jarvis played in the 1971 Quebec International Pee-Wee Hockey Tournament with a minor ice hockey team from Toronto.

Selected by the Washington Capitals in the 1978 NHL Entry Draft, Jarvis also played for the Minnesota North Stars, Los Angeles Kings, and Toronto Maple Leafs. Jarvis was the winner of the 1978–79 Gary F. Longman Memorial Trophy, which is awarded to the player voted to be most outstanding in his first season in the International Hockey League by the league coaches. During the 1982–83 season, Jarvis won the Phil Esposito Trophy, which is awarded to the leading scorer of the Central Hockey League for the regular season. His last season of professional hockey was in 1989–90 with the Newmarket Saints. Jarvis was the head coach of the Newmarket Hurricanes of the OPJHL and was an assistant coach with the Barrie Colts for three seasons.

Personal life
Now living in Ontario with his wife (Darlene) and four daughters (Darcie, Corie, Terrie and Leslie), Wes and former NHL teammate Mike Gartner own and run National Training Rinks, a hockey training and instruction facility. He is the uncle of Alex Foster. He is cousin with Doug Jarvis.

References

External links

Profile at hockeydraftcentral.com

Category:1958 births
Category:Canadian ice hockey centres
Category:Hershey Bears players
Category:Living people
Category:Los Angeles Kings players
Category:Minnesota North Stars players
Category:Sportspeople from Toronto
Category:Sudbury Wolves players
Category:Toronto Maple Leafs players
Category:Washington Capitals draft picks
Category:Washington Capitals players
Category:Windsor Spitfires players


Furniture Wheels

To make the furniture more and more beautiful, stylish and decorative, the furniture wheels are made with equal and similar color combinations. The furniture wheels are specific part of the furniture casters, which are made with furniture wheels, fitted in them. The frames of the these furniture casters are made with the comparatively strong, hard, and tough material, such as, the blue or black steel, aluminum or the hard plastic. But, sometimes, the furniture wheels are made with the different type of the material, such as, PVC plastic, the fiber glass, the synthetic or natural rubber and sometimes with bras.

The different material to make the wheels makes the furniture casters more and more attractive to see. The PVC, fiberglass, synthetic or natural rubber made furniture wheels UK made in different colors to add beauty to the wheel casters and the furniture as well. The bras wheel casters and the furniture wheels are used in the high quality and luxury furniture as they may be polished again and again to give sparkling shine.

Various types and kinds of the furniture wheels UK are in use for the different types of the furniture. For the heavy duty furniture wheels, the heavy duty wheel casters are used. The heavy duty furniture wheels for the heavy duty casters are also made with the material of similar quality and durability. Normally, the heavy duty furniture wheels are made for the racks and cupboards for the industrial and commercial usage. They are made with very strong and hard blue or black steel or aluminum in order to sustain the weight of even tones. The heavy duty furniture wheels for the industrial or commercial usage are made with steel or other hard metal as well. They are protected from the environmental rusting or the moisture with the coating of galvanization or any other type of electroplating treatment.


Staxel Reporter – Issue #30

Hey everybody, in this issue of the Staxel Reporter we’ll be talking about biomes and take a look at the new pets. Enjoy!

Biomes

Many of you will be excited to hear that we’ve recently begun work on adding different biomes to the Staxel world! Let’s take a look at what’s in the works.

The forest biome is the first biome we started working on, and as we only recently started working on biomes it’s currently the only one that’s screenshot ready. Right now it doesn’t offer anything other than a pretty zone on your island, but we aim to give different biomes real purpose so in the future you might find cool stuff in there that you won’t find anywhere else!

The next one in the works is likely going to be more of a wide open space, something along the lines of a meadow or a savannah. We’re very excited for the new biomes as it will be a great way to break up the island. Now, we need to do something about the ocean of flowers…

We’ll definitely be posting pics when the next biome is ready so make sure you keep an eye on any of the following:

New Pets

The variety of pet breeds in Staxel has just grown! Before, at the end of the tutorial you could only recieve one type of breed for both cats and dogs. Now, however, you’re rewarded with a pet of a random breed! But who cares about what I’m saying? Let’s look at some cute pet pictures!

D’awww, just look at these adorable dogs! The new breeds are bloodhounds, dalmation, and border collies. There may be more breeds in the future but these should do for now at least! Now let’s put the dogs safely away and bring out the cats.

You have no idea how hard it was to get them to line up for this picture… does anyone have any antiseptic? Anyway, here they are! We’re hoping these extra breeds should help your pet feel a bit more unique to you.

That about sums it up for the new pets. We’ll certainly be adding new pets in the future other than cats and dogs, so there’s more to look forward to in the pet department!

That’s it for this issue of the Staxel Reporter. You can grab ‘Staxel – Sprout Edition’ at www.playstaxel.com/store.

Thanks for reading!

–EGadd


Q:

calculus continuity function of unknown

I am watching a video on Introduction of Continuity, and I am confused by this statement made by the professor: "when $x$ is lesser than $1$, it wants to be $ax - 4$, but I want, when $x$ is $1$, it to be $a\cdot 1 - 4 = 3$, which is equal to $a - 4 = 3$, which implies $a$ is $7$."
Since $x$ is $1$, is satisfies the condition $x \geq 1$, and it should evaluate to $3(1)^2$ instead of $a(1)-4=3$.
Can anyone explain?

A:

The problem we're trying to solve is:

Find all values of $a$ such that $f$ is continuous everywhere.

In particular, we need $f$ to be continuous at $x = 1$. By definition, this means that:

$$
\lim_{x \to 1} f(x) = f(1) = 3 \cdot 1^2 = 3
$$

For the two-sided limit to equal some number (in this case, $3$), we at least need the two-sided limit to exist. To guarantee existence, each of the one-sided limits must also exist and match each other. So in particular, we know that:

$$
\lim_{x \to 1^-} f(x) = 3
$$

But if $x$ is approaching $1$ from the left, then $x < 1$, so we know which piece of $f$ to use. This means that:

$$
\lim_{x \to 1^-} (ax - 4) = 3 \iff a \cdot 1 - 4 = 3 \iff a = 7
$$


Q:

HMAC-SHA1 in Rust

I'm trying to apply HMAC-SHA1 in order to check some content but I'm unable to make it work.
These are the tests I have:
#[cfg(test)]
mod tests {

    use crypto::hmac::Hmac;
    use crypto::mac::Mac;

    use crypto::sha1::Sha1;
    use std::str::from_utf8;

    const BODY_CONTENT: &'static str = r#"bodystring"#;
    const KEY: &[u8] = b"secret_key";
    const COMPUTED_HMAC: &'static str = "97049623b0e5d20bf6beb5313d80600e3d6abe56";

    #[test]
    fn test_hmac_sha1() {
        let mut mac= Hmac::new(Sha1::new(), KEY);
        mac.input(BODY_CONTENT.as_bytes());
        let result = mac.result();
        let code = result.code();
        assert_eq!(COMPUTED_HMAC.as_bytes(), code);
        assert_eq!(COMPUTED_HMAC, from_utf8(&code).unwrap_or("failed"));
    }

    #[test]
    fn test_hmac_sha1_direct() {
        let hash = hmacsha1::hmac_sha1(KEY, BODY_CONTENT.as_bytes());
        assert_eq!(COMPUTED_HMAC.as_bytes(), hash);
        assert_eq!(COMPUTED_HMAC, from_utf8(&hash).unwrap_or("failed"));
    }
}

I've used this website in order to get the COMPUTED_HMAC by using one string (BODY_CONTENT) and a secret key (KEY).
As you can see, I'm trying to leverage both rust-crypto and hmac-sha1 crates and I obtain the same result with both of them.
The thing is that this result doesn't match with what I get in the website (97049623b0e5d20bf6beb5313d80600e3d6abe56) and the tests fail. You may think that the website is wrong but that's not the case as I'm using it to validate some other hashes generated by Github (I'm working in a Github App) and it works.
Then, obviously, I'm missing some step here but I'm unable to figure it out and I would really appreciate your help.

A:

The correct hash is returned, it's just not in the representation you expected. The hash is returned as raw bytes, not as bytes converted to ASCII hexadecimal digits.
If we print the hash code array as hex, like this:
println!("{:02x?}", code);

then we can see that it matches your string:
[97, 04, 96, 23, b0, e5, d2, 0b, f6, be, b5, 31, 3d, 80, 60, 0e, 3d, 6a, be, 56]
// 97049623b0e5d20bf6beb5313d80600e3d6abe56

whereas the string "97049623b0e5d20bf6beb5313d80600e3d6abe56" looks like this:
[39, 37, 30, 34, 39, 36, 32, 33, 62, 30, 65, 35, 64, 32, 30, 62, 66, 36, 62, 65,
 62, 35, 33, 31, 33, 64, 38, 30, 36, 30, 30, 65, 33, 64, 36, 61, 62, 65, 35, 36]

Using itertools, we can convert the former to the latter like this:
assert_eq!(
    COMPUTED_HMAC,
    code.iter().format_with("", |byte, f| f(&format_args!("{:02x}", byte))).to_string());