function [w, b] = sepcore( in, w, b, ep, p, h, mth)

% set default method to bell
if nargin == 6
	  mth = 0;
end

% remember and setup some stuff
mm = floor(length(in)/p);

% count time
tic

% Bell & Sejnowski method
%
%   Working command: bsep( smix, 2, 20, .05, 0);
%
%       Comments:
%               As decribed in: "An Information-Maximization Approach
%               to Blind Separation and Blind Deconvolution", 
%                       from Neural Computation, 7, 6 Nov(?) 95.
%       Works good.  I'm not sure whether I should keep the biases ...
%       But it works great!
%
%       The commented out delta rule is Bell's original.
%       The one right under it is Bell's rule after Amari's suggestions
im = eye( rows( in));
if mth == 0
        on = ones( 1, p);
        for i = 1:ep
                for k = 1:p:mm*p
                        x = in( :, k:k+p-1);
                        u = w * x;
                        y = ptanh( u + b' * on);
                        %dw = p * inv( w') - (2 * y * x');
                        dw = p * ( im - 2 * y * u') * w;
                        db = sum( -2 * y');
                        w = w + h * dw;
                        b = b + h * db;
                end
        end

% Amari, Cichocki & Young method
%
%   Working command: bsep( smix, 2, 20, .05, 1);
%
%   Comments:  
%               As described in "A New Learning Algorithm for
%               Blind Signal Separation", from NIPS 8 1996.
%       Works pretty good.  Vector fy can be any monotonic function.
%       The big and weird one (commented out) is in the paper 
%       but it sucks. Can also use the identity function.  The tanh function
%       is pretty good too but twice as slow.

elseif mth == 1
        im = eye( rows( in));
        for i = 1:ep
                for k = 1:p:mm*p
                        y = w * in( :, k:k+p-1);
%                       fy = .75*(y.^11)+6.23*(y.^9)-4.6667*(y.^7)...
%                                       -11.75*(y.^5)+14.5*(y.^3);
                        fy = ptanh( y);
                        dw = (im - fy * y') * w;
                        w = w + h * dw;
                end
        end

% Cichocki, Kasprzak & Amari method
%
%       Working command:  bsep( 100*smix, 2, 20, .000000001, 2);
%
%       Comments:
%               As described in: "Local Adaptive Learning Algorithms
%               for Blind Separation of Natural Images" from ??(I forget right now...)
%
%       Right now I'm using the global learning rule.  Maybe I'll
%       implement the local one as well ...

elseif mth == 2
        im = eye( rows( in));
        y = in;
        for i = 1:ep
                for k = 1:p:mm*p
                        y = w * in( :, k:k+p-1);
                        gy = y.^3;
                        fy = y;
                        dw = (im - fy * gy') * w';
                        w = w + h * dw;
                end
        end
end

% report time
toc
