function [varargout]=signsubsets(phi,varargin)
%
% This function designs a weight function for inner products of the form:
%
% <f,g> = integrate(f*g*w);
%
% The weight funtion w is constructed such that the absolute value of the
% inner product of f and and g is minimized, so that w remains positive.
% The MATLAB function lsqnonneg (least-squares & non-negative) is
% implemented to do this.
% The user has the option of creating a basis orthogonal to g whose first
% element is f.
%
% INPUT
%     f: an Nx1 vector, taken to be the data.
%     g: an Nx1 vector, same size as signal.
%------------------------------------------------------------------------
% NOTE: The product of signal and artifact must have the same sign over a
% subset fo values, and differ in sign over a subset of values.
%------------------------------------------------------------------------
%     t: (optional) variable of integration in inner product.
% basis: (optional) a matrix of linearly indepenet waveforms used to
%        create a basis.  This basis is constructed out of the subspace
%        of functions orthogonal to user input 'g'.
%
%
% Also available as part of VSPLAB from site:
% http://www.ess.washington.edu/~joshuadc

if(nargin==0)
    varargin{1}=0;
    varargin{2}=0;
end;
if(nargin==1);
    varargin{2}=0;
end;

[Vmax,Imax]=max(phi);
[Vmin,Imin]=min(phi);
%give the weight function the same signal width as a phi
Isame=find(phi>=0);
Ip=Imax;
Ipcnt=zeros(length(phi),1);
k=0;
while(~isempty(intersect(Ip,Isame)))
    k=k+1;
    Ip=Imax-k;
    Ipcnt(k)=Ip;
end;
k=0;
In=Imax;
Incnt=zeros(length(phi),1);
while(~isempty(intersect(In,Isame)))
    k=k+1;
    In=Imax+k;
    Incnt(k)=In;
end;
Ipint=union(Ipcnt,Incnt); %get continuous subset over which phi>0
Ipint=union(Imax,Ipint); %include the maximum value index
Ipint=setdiff(Ipint,0); %get rid of zero
Ipint=setdiff(Ipint,min(Ipint)); %get rid of initial false value
Ipint=setdiff(Ipint,max(Ipint)); %get rid of last false value
%FIND SUBSET OF CONTINUOUSLY NEGATIVE VALUES
Idiff=find(phi<0);
Ip=Imin;
Ipcnt=zeros(length(phi),1);
k=0;
while(~isempty(intersect(Ip,Idiff)))
    k=k+1;
    Ip=Imin-k;
    Ipcnt(k)=Ip;
end;
k=0;
In=Imin;
Incnt=zeros(length(phi),1);
while(~isempty(intersect(In,Idiff)))
    k=k+1;
    In=Imin+k;
    Incnt(k)=In;
end;
Inint=union(Ipcnt,Incnt); %get continuous subset over which phi>0
Inint=union(Imin,Inint); %include the maximum value index
Inint=setdiff(Inint,0); %get rid of zero
Inint=setdiff(Inint,min(Inint)); %get rid of initial false value
Inint=setdiff(Inint,max(Inint)); %get rid of last false value

if(nargout==1);
    varargout{1}=Ipint;
end

if(nargout==2);
    varargout{1}=Ipint;
    varargout{2}=Inint;
end

if(nargout==0||(strcmp(varargin{1},'plot')||strcmp(varargin{2},'plot')))
    plot(phi,'linewidth',1.6); hold on;
    s=plot(Ipint,zeros(size(Ipint)),'r'); set(s,'linewidth',2);
    s=plot(Inint,zeros(size(Inint)),'k'); set(s,'linewidth',2);
    s=title('Function with Continous Subsets that Include Global Extrema');
    set(s,'fontsize',14);
    set(gca,'fontsize',14);
    legend('Input','Maximum','Minimum');
    set(gcf,'Color', [1,1,1]);
    set(gca,'Color', [0.85,0.9,0.95]);
end;