% NSCI850_solveRBFnet.m
%
% computes ,ulti-sensory integration and reference frame transformation
% within the same distibuted code
%
% created by G. Blohm for NSCI 850 (2018)

clear;
warning off;

%% creating training set
%-----------------------

% defining the exact input-output relationship
N = 10000; % length of training set
ff = 100; % firing rate (for scaling)
input(1,:) = 15*(randn(1,N)); % random set of retinal hand position (in visual degrees)
input(2,:) = input(1,:) + 10*(randn(1,N)); % random set of associated proprioceptive hand position
VarVis = (abs(1*(randn(1,N)+2))+1).^2; % visual variance
VarPro = (abs(1*(randn(1,N)+3))+1).^2; % proprioceptive variance
input(3,:) = input(2,:) + 10*(randn(1,N)); % eye position
VarEye = (abs(1*(randn(1,N)+3))+1).^2; % eye position variance
VarTot = (VarVis+VarEye).*VarPro./(VarVis+VarEye+VarPro); % total variance (Bayesian)
output = VarTot.*((input(1,:)+input(2,:))./(VarVis+VarEye) + input(2,:)./VarPro); % optimal spatial target position (Bayesian)

%% initialize network parameters (RFs, gain fields)
Nh = 100; % number of neurons
% below are random initializations of RF centers & width as well as gain
% modulation properties of each neuron
xx = 200*(rand(Nh,1)-0.5);
ee = 2*(rand(Nh,1)-0.5);
pp = 2*(rand(Nh,1)-0.5);
ees = 10*(rand(Nh,1)-0.5);
pps = 10*(rand(Nh,1)-0.5);
xxs = 10*(rand(Nh,1)-0.5);
oo = 1*(rand(Nh,1)-0.5);

%% compute read-out weights for RBF net
% variance modulates gain field strength inversely as in Ma, Beck & Pouget
% (2006)
g_ep = repmat(oo,1,length(input(3,:))) + ee*(input(3,:))+pp*(input(2,:))+xxs*(1./VarVis)+ ees*(1./VarEye)+pps*(1./VarPro); % gain field
RF_x = exp(-(repmat(xx,1,length(input(3,:)))-repmat(input(1,:),Nh,1)).^2./(2*100)); % receptive field
aact = g_ep.*RF_x; % gain flield * RF

w=pinv(aact')*output'; % inverse computation of read-out weights
output_test = aact'*w;

%% plot network output vs. desired output
figure; 
plot(output,output_test,'.')
xlabel('Desired output (deg)')
ylabel('RBFnet output (deg)')

