clear all
close all
baseline = 5;
noisiness = 4;
oomph = 1.5;
%1 Understanding least squares regression and gradient descent 

%Conceptual understanding:
%Say we have an experiment with 3 trials - this is nice because
%two would be trivial and we can visualize 3. 

xLuminance = [1:3]'; %Arbitrary units. Column vector of IV
ySpikes = oomph * xLuminance + baseline + randn(length(xLuminance),1).*noisiness; %Noisy integrate and fire

%Plotting this to see what is going on
figure
h = plot(xLuminance,ySpikes,'o') %Put a handle on things
h.MarkerSize = 18;
h.MarkerFaceColor = 'b';
h.MarkerEdgeColor = 'k';
xlabel('Luminance (arbitrary units)')
ylabel('Spikes in ips')
set(gca,'Fontsize',18)
xlim([0.5 3.5])

%% Looks legit, but we need to now implement the equation from the lecture
%We need to interpret the experiment visually
%Each trial will now be a dimension in the experiment-space
%Plot the entire experiment (both IVs and DVs) as vectors
figure
h1 = plot3([0, xLuminance(1)], [0, xLuminance(2)], [0, xLuminance(3)])
h1.Color = 'b';
h1.LineWidth = 2; 
hold on
h2 = plot3([0, ySpikes(1)], [0, ySpikes(2)], [0, ySpikes(3)])
h2.Color = 'g';
h2.LineWidth = 2;
legend('Luminance','Spikes')
xlabel('trial 1')
ylabel('trial 2')
zlabel('trial 3')
set(gca,'fontsize',18)
rotate3d on

%View on the same data in 3d - now we can use linear algebra. Instead of
%having 3 points in 2d, we have 2 points in 3d. 

%Now let's find the scalar, that - when multiplied with the IV vector,
%minimizes the distance to the measurements. This scalar is called beta.

beta = dot(ySpikes,xLuminance)/dot(xLuminance,xLuminance); %Normalized projection
prediction = beta*xLuminance; %Simplest possible case

%Add this to the plot - to thicken it
h3 = plot3 ([0, prediction(1)], [0 prediction(2)], [0 prediction(3)]);
h3.LineWidth = 2;
h3.Color = 'k';
h3.LineStyle = ':';
legend('Luminance','Spikes','Prediction')

%Let's explicitly add the distance between the error
h4 = plot3([ySpikes(1), prediction(1)],[ySpikes(2), prediction(2)], [ySpikes(3), prediction(3)]);
h4.LineWidth = 2;
h4.Color = 'r';
legend('Luminance','Spikes','Prediction','Error')
axis equal

%% As per Charlie's suggestion, let's go back to the original plot and add the line thusly defined
%Also, vary beta to see how that changes things
%Let's say we do 20 measurements
maxLuminance = 20;
xLuminance = [0:maxLuminance]'; %A column vector of luminance values
ySpikes = oomph * xLuminance + baseline + randn(length(xLuminance),1).*noisiness; %Noisy integrate and fire
figure
h5 = plot(xLuminance,ySpikes,'o')
h5.MarkerSize = 12;
h5.MarkerFaceColor = 'b';
h5.MarkerEdgeColor = 'k';

%Now that we convinced ourselves that this logic is legit, we can repurpose
%it without visualizing it
beta = dot(ySpikes,xLuminance)/dot(xLuminance,xLuminance); %Normalized projection
prediction = beta*xLuminance; %Simplest possible case

%Let's add the line
pause 
hold on
regressionLineX = linspace(0,maxLuminance,length(xLuminance)); %Equally spaced numbers between 0 and maxLuminance
regressionLineY = beta .* regressionLineX; %That's it
h6 = plot(regressionLineX,regressionLineY,'k')
h6.LineWidth = 2; 

%Add the error
pause
lh = line([xLuminance'; xLuminance'],[prediction'; ySpikes'], 'color','r')
pause
% Now: Multiple regression - multiple predictors
%Each predictor will be a column in the matrix
%Simplest case: Adding an offset. Even a constant will turn this into a
%multiple regression problem 
%To represent the baseline, we need to add a constant value per trial
offset = ones(length(xLuminance),1); %Representing the baseline firing rate in absence of visual stimulation
designMatrix = [xLuminance, offset]; 
%The equation calls for an inverse, we will use the pseudoinverse to do the
%best we can
betas = pinv(designMatrix' * designMatrix) * designMatrix'*ySpikes; %This should yield 2 betas, one updated slope and an intercept

prediction = designMatrix * betas; %Same as before, but X is no longer a vector, it is a matrix

%Let's just plot the new regression line on top of the old one
%regressionLineYBar = betas(1) * designMatrix(:,1) + betas(2) ...
%    * designMatrix(:,2); 
plot(regressionLineX,prediction,'color','m','linewidth',2)
%pause
%plot(regressionLineX,prediction,'color','y','linewidth',2) %We could also just use the prediction we made

%% 2 Dropping marbles - which metric to find beta is the right one? 
%Let's explore the space of betas
startExploration = beta-2;
endExploration = beta+2;
numBeta = 200; %Fine resolution
testBetas = linspace(startExploration, endExploration, numBeta);
for ii = 1:numBeta %Go through all the betas
    prediction = testBetas(ii) .* xLuminance;
    %We need to now explicitly represent the distance metric
    errorMetric(ii) = sum((prediction-ySpikes).^2); %Sum of squares
end

%Plot that
figure
plot(testBetas,errorMetric,'b','linewidth',3)
line([beta beta], [min(ylim) max(ylim)], 'color','r');
xlabel('Beta')
ylabel('Error metric')
set(gca,'fontsize',18)
set(gcf,'color','w')


%% 3 Going deeper into optimization - and gradient descent
beta
startExploration = beta-2;
endExploration = beta+2;
testBetas = linspace(startExploration,endExploration,200);
for ii = 1:length(testBetas)
    prediction(:,ii) = testBetas(ii)*xLuminance;
    distanceSum(1,ii) = sum(prediction(:,ii)-ySpikes); %Simple
    distanceSum(2,ii) = sum((prediction(:,ii)-ySpikes).^2); %Squared
    distanceSum(3,ii) = sum(abs(prediction(:,ii)-ySpikes)); %Absolute value
    distanceSum(4,ii) = sum(log(1 + (prediction(:,ii)-ySpikes).^2)); %Lorentzian
    distanceSum(5,ii) = sum(log(prediction(:,ii)-ySpikes)); %Weirdness
    distanceSum(6,ii) = sum(log(prediction(:,ii)-ySpikes).^2); %Weirdness

end

figure
plot(testBetas,distanceSum(2,:),'b','linewidth',3)
line([beta beta],[min(ylim) max(ylim)],'color', 'r')
xlabel('Beta')
ylabel('Distances')
set(gca,'fontsize',18)

figure
for ii = 1:6
    subplot(2,3,ii)
plot(testBetas,distanceSum(ii,:),'b','linewidth',3)
line([beta beta],[min(ylim) max(ylim)],'color', 'r')
xlabel('Beta')
ylabel('Distances')
if ii == 1
    title('Summed differences')
elseif ii == 2
    title('Sum of squared differences')
elseif ii == 3
    title('Sum of absolute differences')
elseif ii == 4
    title('Lorentzian?')    
elseif ii == 5
    title('Sum of log differences')        
elseif ii == 6
    title('Sum of squared log differences')    

end
    
set(gca,'fontsize',18) 

end
set(gcf,'color','w')
%% 4 Gradient descent revisited
figure
for ii = 1:length(testBetas)
    subplot(1,2,1)
    hold on
    prediction = testBetas(ii)*xLuminance;
    h7 = plot(xLuminance, ySpikes, 'o');
    set(h7,'markersize',12)
    set(h7,'markerfacecolor','b')
    set(h7,'markeredgecolor','k')
    xlabel('Luminance in cd/m^2')
    ylabel('Spikes in ips')
    set(gca,'fontsize',18)
    regressionLineX = linspace(0,20,100); %Intrapolate

    regressionLineY = testBetas(ii) * regressionLineX;
    h8 = plot(regressionLineX, regressionLineY, 'k');
    set(h8,'linewidth',2);
    
for jj = 1:length(ySpikes)
  h(jj) = line([xLuminance(jj) xLuminance(jj)],[prediction(jj) ySpikes(jj)],'color','r');
end
    ylim([-5 1.5*max(ySpikes)])
    subplot(1,2,2)
    hold on
    plot(testBetas(ii),distanceSum(2,ii),'.','color','b','linewidth',3)
    xlabel('Beta')
    ylabel('Distance')
    set(gca,'fontsize',18)
    xlim([startExploration endExploration])
    ylim([0 max(distanceSum(2,:))])
    shg
    pause(0.1)
    if ii < length(testBetas)
    delete(h8)
    delete(h)
    end
end


%% 5 PCA
corr = 0.9 %Possible values: 0, 0.7, 0.9, 0.3, -0.7 
A = randn(1000,2); %uncorrelated data

if corr == 0.7 % Corr ~ 0.7
B = randn(1000,1); 
C = randn(1000,1);
D = (B + C)./2; %correlated data (B and D)
elseif corr == 0
B = randn(1000,1); 
D = randn(1000,1);
elseif corr == 0.9
B = randn(1000,1); 
C = randn(1000,1);
D = (3.*B + C)./4; %correlated data (B and D)
elseif corr == 0.3
Q = randn(1000,1);     
B = randn(1000,1); 
C = randn(1000,1);
D = ((3*Q) + B + C)./5; %correlated data (B and D)
elseif corr == -0.7
B = randn(1000,1); 
C = randn(1000,1);
D = (B + C)./2; %correlated data (B and D)
D = D.*-1;
end

dotVec = [B,D];
dotProducts = zeros(1000,1);
spreads = zeros(360,1);

fullScreen
subplot(1,2,1)
plot(A(:,1),A(:,2),'.')
xlim([-5 5])
ylim([-5 5])

fullScreen
subplot(1,3,1)
plot(B(:,1),D(:,1),'.')
xlabel('x Value')
ylabel('y Value')
temp = corrcoef(B,D); 
title(['corr = ', num2str(temp(1,2))])
xlim([-5 5])
ylim([-5 5])
M = [5;0];
hold on
plotv(M,':');
plotv(-1.*M,':');
h = findobj('linestyle',':');
set(h,'color','r')
set(h,'linestyle','-')
for ii = 1:1000
    dotProducts(ii,1) = dot(dotVec(ii,:),M);
end

subplot(1,3,2)
histogram(dotProducts,25)
xlabel('Dot product')
ylabel('#')
xlim([-20 20])
ylim([0 150])


subplot(1,3,3)
plot(1,std(dotProducts),'.','color','k');
xlabel('Vector angle')
ylabel('Variation')
xlim([0 361])
ylim([0 7])
set(gca,'XTick', [0 200 360])
hold on
%set(h,'color','r')
%columnVector = [-5 5]'
%hold on

pause
%R = [0 -1; 1 0];
for aa = 1:360
subplot(1,3,1)
angle = deg2rad(aa);
R = [cos(angle) -sin(angle); sin(angle), cos(angle)];
RM = R*M; 
set(h,'linestyle','none')
plotv(RM,':');
plotv(-1.*RM,':');
h = findobj('linestyle',':');
set(h,'color','r')
set(h,'linestyle','-')
pause(0.005)

for ii = 1:1000
    dotProducts(ii,1) = dot(dotVec(ii,:),RM);
end

subplot(1,3,2)
histogram(dotProducts,25)
xlabel('Dot product')
ylabel('#')
xlim([-20 20])
ylim([0 150])

subplot(1,3,3)
spreads(aa) = std(dotProducts);
plot(aa,spreads(aa),'.', 'color','k');
xlabel('Vector angle')
ylabel('Variation')
set(gca,'XTick', [0 200 360])

hold on
%columnVector2 = R*columnVector;
%plot(columnVector2,'color','r')
%myVector = [-5 0; 5 0]l
%myVector2 = R*myVector;
%line(myVector(:,1),myVector(:,2),'color','r')
%pause
%line(myVector2(:,1),myVector2(:,2),'color','b')
end
maxSpreads = max(spreads);
temp = find(spreads == maxSpreads);
title(['Max angle = ', num2str(temp(1))])
line([temp(1) temp(1)],ylim,'color','r')

subplot(1,3,1)
M = [5;0];
angle = deg2rad(temp(1,1)-1);
R = [cos(angle) -sin(angle); sin(angle), cos(angle)];
RM = R*M; 
set(h,'linestyle','none')
plotv(RM,':');
plotv(-1.*RM,':');
h = findobj('linestyle',':');
set(h,'color','g')
set(h,'linestyle','-')
set(h,'linewidth',3)

%% 6 Introducing while statements and mixing if statements into loops
%Let's find 10 orthogonal vectors by random trial and error
orthoCounter = 1 ; %This is our counter of orthogonal vectors
loopCounter = 0; %This counts the number of loops it took to get them
while orthoCounter < 11 % We don't know how many tries we need, so while is suitable
    A = randn(1,2); %Draw a 2d row vector randomly from a normal distribution
    B = randn(2,1); %Draw a 2d column vector randomly from a normal distribution
    testThem = A*B; %Take the inner product
    loopCounter = loopCounter + 1; 
    if abs(testThem) < 0.000001 %Ideally, we would IPs that are exactly 0, due to MATLABs numerical precision, this will never happen, so we need to check for close enough to zero
        disp('Orthogonal!') %Expression of joy
        rowContainer(orthoCounter,:) = A; %Capture the row vector
        columnContainer (:,orthoCounter) = B; %Capture the corresponding column vector. We should initialize these
        orthoCounter = orthoCounter + 1; %Increment the index
    end %Close the if statement #1
    
end %Close the while loop
%The reason I'm showing the next command is to give an example of how to
%mix characters and numbers in the same statement.
disp(['It took ', num2str(loopCounter), ' repeats to find 10 orthogonal vectors'])

%% 7 Now that we found them, let's plot them
figure %Unless you call this explicitly and there is not one open, it will open one for you. Here, it is good to open one explictly
%If you have more than one figure you draw to, it helps to open new figures
for ii = 1:10
    plot([0 rowContainer(ii,1)],[0 rowContainer(ii,2)],'color','b') %Plot x-piece of row vector from origin, then y-piece in blue
    hold on %Put hold on, otherwise the 2nd vector will overwrite the first one
    plot([0 columnContainer(1,ii)], [0 columnContainer(2,ii)],'color','r') %Plot x-piece of column vector from origin, in red
    hold off %If I don't put the hold off again, the next iteration of the loop will all draw them on top of each other
    %Exercise: Play around with hold to see what will happen
        axis square
        axis equal %If axis units are not equal, it will look distorted
        pause %Only move on to the next loop on button press

end

%% 8 Putting it *all* together - matrix multiplication, loops, random numbers
%We want to create a constellation of dots that represents some data or
%system
%We want to rotate it with a orthogonal matrix of unit vectors. These
%rotate the vectors without changing their lengths 
%close all %This closes all pre-existing figures
figure %Open a new figure
A = randn(1000,2); %Drawing 1000 random vectors in 2 dimensions from a normal distribution
theta = -1 / 180 * pi; %1 degree in radians
B = [cos(theta) sin(theta); -sin(theta) cos(theta)]; %One possible orthogonal unit vector
for ii = 1:359
    plot(A(:,1),A(:,2),'.','markersize',20); %Plot the endpoints of the vectors, not the vectors themselves. As dots
    xlim([-4 4]) %This will capture all dots from a normal distribution in all likelihood
    ylim([-4 4])
    pause(0.01) %With an argument, pause tells MATLAB how long to wait before moving on, in seconds
    
    A = A*B; %Updating A each iteration of the loop with the rotation from the ortogonal matrix
    shg
end


