clear 
clc

%single intervention
%generate a random A matrix with n_s sectors n_countries
n_s=200;
n_c=49;
%random matrix generation or load exiobase A matrix
%load the A,b, and y for EXIOBASE
load('A.mat')
load('b.mat')
load('y.mat')
% c = 0;
% b = 0.001; %max of the tech coef
% A = (b-c).*rand(n_s*n_c,n_s*n_c) + c;
% 
% b=rand(1, n_s*n_c);
% y=rand(n_s*n_c,1 );
L= (eye(n_s*n_c,n_s*n_c)-A);
%Assume you know the initial Leontief inverse
inv_L=inv(L); 


 %multiple interventions: alpha create the table of interventions number of
 %rows =alpha, number of columns =3 corresponding to the values of
 %k,i,  c_i, 
 k=150; %the column number of the construction work
 i=[104;108;150]; %the input of k that are subjected to change
 c_i=[-0.4;-0.45;0.11]; %the amount of the change in the inputs
 interventions_table=zeros(n_c*n_c,size(c_i,1));
 for t=1:size(c_i,1)
  inputs=k+0:200:9800;
  outputs=i(t)+0:200:9800;
 interventions_table( 1+n_c*n_c*(t-1) :n_c*n_c +n_c*n_c*(t-1) ,1:2) = (allcomb(inputs,outputs));
 interventions_table( 1+n_c*n_c*(t-1) :n_c*n_c +n_c*n_c*(t-1) ,3)=repmat(c_i(t,1),n_c*n_c,1);
 end
   
 
 interventions_set=( unique(interventions_table(:,1)));
  alpha=size( unique(interventions_table(:,1)),1); %number of interventions in different sectors
  n=size(A,1); 
  u=zeros(n,alpha);
  v=zeros(alpha,n);
  i=eye(alpha,alpha);
  inv_L=inv(eye(n,n)-A);
  y=transpose(y);
 
  %calculation with woodbury
  for t=1:alpha
  v(t,interventions_set(t,1))=1;
  dum=interventions_table(:,1)==interventions_set(t,1);
  u(interventions_table(dum,2),t)=A((interventions_table(dum,2)),interventions_set(t,1)).*interventions_table(dum,3);
  end
 
  tic
  delta_r_W= b *  inv_L*u*inv(i-v*inv_L*u)*v*inv_L * y;
  toc
     
  %default calculation
  for t=1:size(interventions_table,1)
  A(interventions_table(t,2),interventions_table(t,1))=A(interventions_table(t,2),interventions_table(t,1))+A(interventions_table(t,2),interventions_table(t,1))*interventions_table(t,3);
  end
 
  tic
  new_inv_L=inv(eye(n)-A);
  delta_r_d= b *  (new_inv_L*y - inv_L*y);
  toc
  
  %backslash calculation
  tic
  L= eye(n)-A;
  delta_r_b= b*(L\y - inv_L*y);
  toc
  
  
  
  