sigma = 1;      %sigma determines both the breadth and spacing of the BFs
gabor_sig = 1;  %sigma of gaussian in gabor x
gabor_freq = 2; %freq of cosine in gabor x
cos_freq = 3;   %freq of cosine in y
alpha = 0.01;   %learning rate
num_rep=20;  %number of repetitions


x=(-10:0.1:10)';   %input space:  from -10 to 10 by 0.1
c=-10:sigma:10;     %centers of basis functions in each dimension
num_g=length(c);    %number of basis functions in each dimension
g=zeros(length(x),length(c));   %initialization of basis functions

z=exp((-x.^2)/(2*(gabor_sig^2))) .* cos(gabor_freq*x);
%truth:  eval of gabor

for ii=1:num_g,
    for jj=1:length(x),
        g(jj,ii)=exp(-((x(jj)-c(ii))^2/(2*sigma^2)));
        
        %populating evaluation of basis fns
        %NB:  these do not change with learning
    end
end

w=pinv(g)*z;        %pseudoinvert g to find best linear w
z_est=g*w;          %evaluate estimate of y
figure
%then plot
plot(x,z,'r')
hold on
plot(x,z_est,'b')
%%%%%%%%%%%%%%
%%%%


%iterative

w1=0.01*randn(1,num_g);     %start with small random weights

for jj=1:num_rep,
    x1=20*rand-10;          %choose an input at random
    g1=zeros(1,num_g);      %define g1, g evaluated just ounce
    z1=exp( (-x1.^2)/(2*(gabor_sig^2))) .* cos(gabor_freq*x1);
    for ii=1:num_g,
        g1(ii)=exp(-((x1-c(ii))^2)/(2*sigma^2));
        %evaluate functions at this input
    end
    z_est=g1*w1';           %make estimate
    w1=w1+alpha.*g1*(z1-z_est);     %update weights via delta rule (p.1195)
end

%%NB:  below calls upon the full x, y, and g from the top, but uses the
%% iteratively found w1

figure
plot(x,z,'r')
hold on
plot(x,g*w1','b')
