 function varargout = GUI_brenda(varargin)
% GUI_BRENDA M-file for GUI_brenda.fig
%      GUI_BRENDA, by itself, creates a new GUI_BRENDA or raises the existing
%      singleton*.
%
%      H = GUI_BRENDA returns the handle to a new GUI_BRENDA or the handle to
%      the existing singleton*.
%
%      GUI_BRENDA('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in GUI_BRENDA.M with the given input arguments.
%
%      GUI_BRENDA('Property','Value',...) creates a new GUI_BRENDA or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before GUI_brenda_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to GUI_brenda_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help GUI_brenda

% Last Modified by GUIDE v2.5 18-Apr-2009 17:08:09

% Begin initialization code - DO NOT EDIT

gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @GUI_brenda_OpeningFcn, ...
                   'gui_OutputFcn',  @GUI_brenda_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before GUI_brenda is made visible.
function GUI_brenda_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to GUI_brenda (see VARARGIN)

% Choose default command line output for GUI_brenda
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes GUI_brenda wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = GUI_brenda_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%RESLOC%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% The path of input& output file must be specified ..
% The path of input& output file must be specified ..

clear all;

% hhclass=5
% Input data ..The path should be changed based on the location of input..
HHactual = load('C:\G-LUM Finale\Calibration\Input\HH_base.txt'); % input path should be specified ..
HH_lag = load('C:\G-LUM Finale\Calibration\Input\HH_lag.txt'); % input path should be specified ..
EMPactual = load('C:\G-LUM Finale\Calibration\Input\EMP_base.txt');% input path should be specified ..
EMP_lag = load('C:\G-LUM Finale\Calibration\Input\EMP_lag.txt');% input path should be specified ..

%Dimensions of the HH & EMP matrices
dimZone=size(HHactual,1); 
dimHH=size(HHactual,2); 
dimEMP=size(EMPactual,2); 
hhclass=dimHH;
empclass=dimEMP;
unemp=zeros(1,dimEMP); %Assume unemployment rate is zero.

%Fill in 0.0001 when HH or EMP is 0
for i=1:dimZone
    for j=1:dimHH
        if HHactual(i,j) < 0.0001
            HHactual(i,j)=0.0001;
        end;
    end;
    
    for j=1:dimHH
        if HH_lag(i,j) < 0.0001
            HH_lag(i,j)=0.0001;
        end;
    end;
    
    for k=1:dimEMP
        if EMPactual(i,k) < 0.0001
            EMPactual(i,k)=0.0001;
        end;
    end;
    
    for k=1:dimEMP
        if EMP_lag(i,k) < 0.0001
            EMP_lag(i,k)=0.0001;
        end;
    end;

end;
if hhclass>0     HHactual_1=HHactual(:,1); end
if hhclass>1     HHactual_2=HHactual(:,2); end
if hhclass>2     HHactual_3=HHactual(:,3); end
if hhclass>3     HHactual_4=HHactual(:,4); end
if hhclass>4     HHactual_5=HHactual(:,5); end
if hhclass>5     HHactual_6=HHactual(:,6);end

if hhclass>0     HH_lag_1=HH_lag(:,1); end
if hhclass>1     HH_lag_2=HH_lag(:,2); end
if hhclass>2     HH_lag_3=HH_lag(:,3); end
if hhclass>3     HH_lag_4=HH_lag(:,4); end
if hhclass>4     HH_lag_5=HH_lag(:,5); end
if hhclass>5     HH_lag_6=HH_lag(:,6); end
%Calculate zonal total HH
sumHHzone=sum(HH_lag,2); %1*** by 1

%All necessary land in each zone

landrP_lag = load('C:\G-LUM Finale\Calibration\Input\landr.txt'); %residential ...input path must be specified
landbP_lag = load('C:\G-LUM Finale\Calibration\Input\landb.txt'); %basic..input path must be specified
landeP_lag = load('C:\G-LUM Finale\Calibration\Input\landu.txt'); %unusable land...input path must be specified
landuP_lag = load('C:\G-LUM Finale\Calibration\Input\landd.txt'); %undeveloped ..input path must be specified
landcP_lag = load('C:\G-LUM Finale\Calibration\Input\landc.txt'); %commerical...input path must be specified
landsP_lag = load('C:\G-LUM Finale\Calibration\Input\lands.txt'); %Streets & Highways...input path must be specified

for i=1:dimZone
    if(landrP_lag(i) < 0.00001) landrP_lag(i) = 0.00001;end;
    if(landbP_lag(i) < 0.00001) landbP_lag(i) = 0.00001;end;
    if(landeP_lag(i) < 0.00001) landeP_lag(i) = 0.00001;end;
    if(landuP_lag(i) < 0.00001) landuP_lag(i) = 0.00001;end;
    if(landcP_lag(i) < 0.00001) landcP_lag(i) = 0.00001;end;
    if(landsP_lag(i) < 0.00001) landsP_lag(i) = 0.00001;end;
end;

 % The developable land that has already been developed ...
landdRatio = (landrP_lag + landbP_lag  + landcP_lag + landsP_lag )./( landrP_lag + landbP_lag  + landuP_lag + landcP_lag + landsP_lag );  
landdRatio = 1 - (landuP_lag./( landrP_lag + landbP_lag  + landuP_lag + landcP_lag + landsP_lag ));  
ZoneSize = load('C:\G-LUM Finale\Calibration\Input\ZoneSize.txt');%input path must be specified


% landv is the vacant developable land
landv=landuP_lag; 
%landr is the residential land
landr=landrP_lag;

% c denotes the travel time of the region.
c = load('C:\G-LUM Finale\Calibration\Input\TT.txt'); % Travel Time ..Input file path must be specified

%******************************************************************
%******************************************************************
%Calculate "a" (ratio of HH to EMP)
sumHHtype=sum(HH_lag(:,1:dimHH)); %Sum of HH of each type in all zones
sumEMPtype=sum(EMP_lag(:,1:dimEMP)); %Sum of EMP of each type in all zones
a=(ones(1,dimEMP)./sumEMPtype)'*sumHHtype; %dimEMP x dimHH 
if hhclass>0 a_1=a(:,1); end
if hhclass>1 a_2=a(:,2); end
if hhclass>2 a_3=a(:,3); end 
if hhclass>3 a_4=a(:,4); end
if hhclass>4 a_5=a(:,5); end
if hhclass>5 a_6=a(:,6); end

options = optimset('MaxIter',1000000, 'MaxFunEvals', 1000000,'Tolx', 0.01,'TolFun',0.01);
% SUmala modified the option for Tolx
% House Category 1 Calibration ...

%Initialization ...
% Different intial starting must be tried out to overcome the local optimality problem
%eta=0.80; alpha= 1.00; beta=-0.005;
%qq=0.01; rr=0.5; ss=0.20;
%bnn = [4.00;1.00;-1.00;-2.00];

if hhclass>0
    
iv=load('C:\G-LUM Finale\Calibration\Input\initial_values.txt');
eta=iv(1,1); alpha=iv(1,2); beta=iv(1,3);               %eta=iv(1,1);alpha= iv(1,2);beta=iv(1,3);qq=iv(1,4);rr=iv(1,5);ss=iv(1,6);bnn(1)=iv(1,7).....bnn(4)=iv(1,10)
qq=iv(1,4); rr=iv(1,5); ss=iv(1,6);
bnn = iv(1,7:6+dimHH)';
%'hellp'
x0=vertcat(eta,alpha,beta,qq,rr,ss,bnn);
y = @(x) NLLS_RESLOC(x,HHactual_1,HH_lag_1,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_1,c,dimZone,dimHH,dimEMP,HH_lag);
[x,fval,exitflag]= fminsearch(y,x0,options);
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_1_parameters.dat',x); % output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_1_fval.dat',fval);% output file path must be specified

end

%strresponse= input('reached here','s');

%******************************************************************
%******************************************************************
if hhclass>1
% House Category 2 Calibration ...
%eta=0.80; alpha= 1.00; beta=-0.05;
%qq=0.01; rr=0.5; ss=0.20;
%bnn = [1.00;4.00;1.00;-2.00];
eta=iv(2,1); alpha=iv(2,2); beta=iv(2,3);               %eta=iv(1,1);alpha= iv(1,2);beta=iv(1,3);qq=iv(1,4);rr=iv(1,5);ss=iv(1,6);bnn(1)=iv(1,7).....bnn(4)=iv(1,10)
qq=iv(2,4); rr=iv(2,5); ss=iv(2,6);
bnn = iv(2,7:6+dimHH)';
x0=vertcat(eta,alpha,beta,qq,rr,ss,bnn);
y = @(x) NLLS_RESLOC(x,HHactual_2,HH_lag_2,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_2,c,dimZone,dimHH,dimEMP, HH_lag);
[x,fval,exitflag]= fminsearch(y,x0,options);
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_2_parameters.dat',x);% output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_2_fval.dat',fval);% output file path must be specified
end

% House Category 3 Calibration ...
if hhclass>2 
% eta=0.80; alpha= 1.00; beta=-0.05;
% qq=0.01; rr=0.5; ss=0.20;
% bnn = [-2.00;1.00;4.00;1.00];
eta=iv(3,1); alpha=iv(3,2); beta=iv(3,3);               %eta=iv(1,1);alpha= iv(1,2);beta=iv(1,3);qq=iv(1,4);rr=iv(1,5);ss=iv(1,6);bnn(1)=iv(1,7).....bnn(4)=iv(1,10)
qq=iv(3,4); rr=iv(3,5); ss=iv(3,6);
bnn = iv(3,7:6+dimHH)';
x0=vertcat(eta,alpha,beta,qq,rr,ss,bnn);
y = @(x) NLLS_RESLOC(x,HHactual_3,HH_lag_3,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_3,c,dimZone,dimHH,dimEMP, HH_lag);
[x,fval,exitflag]= fminsearch(y,x0,options);
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_3_parameters.dat',x);% output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_3_fval.dat',fval);% output file path must be specified

end


if hhclass>3
% House Category 4 Calibration ...
%eta=0.80; alpha= 1.00; beta=-0.05;
%qq=0.01; rr=0.5; ss=0.20;
%bnn = [-2.00;-1.00;2.00;4.00];
eta=iv(4,1); alpha=iv(4,2); beta=iv(4,3);               %eta=iv(1,1);alpha= iv(1,2);beta=iv(1,3);qq=iv(1,4);rr=iv(1,5);ss=iv(1,6);bnn(1)=iv(1,7).....bnn(4)=iv(1,10)
qq=iv(4,4); rr=iv(4,5); ss=iv(4,6);
bnn = iv(4,7:6+dimHH)';
x0=vertcat(eta,alpha,beta,qq,rr,ss,bnn);
y = @(x) NLLS_RESLOC(x,HHactual_4,HH_lag_4,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_4,c,dimZone,dimHH,dimEMP, HH_lag);
[x,fval,exitflag]= fminsearch(y,x0,options);

csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_4_parameters.dat',x);% output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_4_fval.dat',fval);% output file path must be specified
end

if hhclass>4
% House Category 5 Calibration ...
% eta=0.80; alpha= 1.00; beta=-0.05;
% qq=0.01; rr=0.5; ss=0.20;
% bnn = [-2.00;-1.00;2.00;4.00];
eta=iv(5,1); alpha=iv(5,2); beta=iv(5,3);               %eta=iv(1,1);alpha= iv(1,2);beta=iv(1,3);qq=iv(1,4);rr=iv(1,5);ss=iv(1,6);bnn(1)=iv(1,7).....bnn(4)=iv(1,10)
qq=iv(5,4); rr=iv(5,5); ss=iv(5,6);
bnn = iv(5,7:6+dimHH)';
x0=vertcat(eta,alpha,beta,qq,rr,ss,bnn);
y = @(x) NLLS_RESLOC(x,HHactual_5,HH_lag_5,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_5,c,dimZone,dimHH,dimEMP, HH_lag);
[x,fval,exitflag]= fminsearch(y,x0,options);
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_5_parameters.dat',x);% output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_5_fval.dat',fval);% output file path must be specified

end

if hhclass>5
% House Category 6 Calibration ...
%eta=0.80; alpha= 1.00; beta=-0.05;
%qq=0.01; rr=0.5; ss=0.20;
%bnn = [-2.00;-1.00;2.00;4.00];
eta=iv(6,1); alpha=iv(6,2); beta=iv(6,3);               %eta=iv(1,1);alpha= iv(1,2);beta=iv(1,3);qq=iv(1,4);rr=iv(1,5);ss=iv(1,6);bnn(1)=iv(1,7).....bnn(4)=iv(1,10)
qq=iv(6,4); rr=iv(6,5); ss=iv(6,6);
bnn = iv(6,7:6+dimHH)';
x0=vertcat(eta,alpha,beta,qq,rr,ss,bnn);
y = @(x) NLLS_RESLOC(x,HHactual_6,HH_lag_6,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_6,c,dimZone,dimHH,dimEMP, HH_lag);
[x,fval,exitflag]= fminsearch(y,x0,options);
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_6_parameters.dat',x);% output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\RESLOC_6_fval.dat',fval);% output file path must be specified

end

if hhclass >0
% Calculation of the residuals ...
% House Category 1 residual computation ...
x0 =load('C:\G-LUM Finale\Calibration\Output\RESLOC_1_parameters.dat');
y =NLLS_RESLOC_SIM(x0,HHactual_1,HH_lag_1,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_1,c,dimZone,dimHH,dimEMP, HH_lag);
eps_HH1 = load('C:\G-LUM Finale\Calibration\Output\epsilon_RESALLOC.dat');
end

if hhclass>1
% House Category 2 residual computation ...
x0 =load('C:\G-LUM Finale\Calibration\Output\RESLOC_2_parameters.dat');
y = NLLS_RESLOC_SIM(x0,HHactual_2,HH_lag_2,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_1,c,dimZone,dimHH,dimEMP, HH_lag);
eps_HH2 = load('C:\G-LUM Finale\Calibration\Output\epsilon_RESALLOC.dat');
end

if hhclass>2
% House Category 3 residual computation ...
x0 =load('C:\G-LUM Finale\Calibration\Output\RESLOC_3_parameters.dat');
y = NLLS_RESLOC_SIM(x0,HHactual_3,HH_lag_3,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_1,c,dimZone,dimHH,dimEMP, HH_lag);
eps_HH3 = load('C:\G-LUM Finale\Calibration\Output\epsilon_RESALLOC.dat');
end

if hhclass>3
% House Category 4 residual computation ...
x0 =load('C:\G-LUM Finale\Calibration\Output\RESLOC_4_parameters.dat');
y = NLLS_RESLOC_SIM(x0,HHactual_4,HH_lag_4,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_1,c,dimZone,dimHH,dimEMP, HH_lag);
eps_HH4 = load('C:\G-LUM Finale\Calibration\Output\epsilon_RESALLOC.dat');
end

if hhclass>4
% House Category 5 residual computation ...
x0 =load('C:\G-LUM Finale\Calibration\Output\RESLOC_5_parameters.dat');
y = NLLS_RESLOC_SIM(x0,HHactual_5,HH_lag_5,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_1,c,dimZone,dimHH,dimEMP, HH_lag);
eps_HH5 = load('C:\G-LUM Finale\Calibration\Output\epsilon_RESALLOC.dat');
end

if hhclass>5
% House Category 6 residual computation ...
x0 =load('C:\G-LUM Finale\Calibration\Output\RESLOC_6_parameters.dat');
y = NLLS_RESLOC_SIM(x0,HHactual_6,HH_lag_6,EMPactual,sumHHzone,unemp,landv,landdRatio,landr,a_1,c,dimZone,dimHH,dimEMP, HH_lag);
eps_HH6 = load('C:\G-LUM Finale\Calibration\Output\epsilon_RESALLOC.dat');
end

if hhclass>0 eps = [eps_HH1']'; end
if hhclass>1 eps = [eps_HH1';eps_HH2']';end 
if hhclass>2 eps = [eps_HH1';eps_HH2';eps_HH3']';end 
if hhclass>3 eps = [eps_HH1';eps_HH2';eps_HH3';eps_HH4']';end 
if hhclass>4 eps = [eps_HH1';eps_HH2';eps_HH3';eps_HH4';eps_HH5']';end 
if hhclass>5 eps = [eps_HH1';eps_HH2';eps_HH3';eps_HH4';eps_HH5';eps_HH6']';end 
% Integrating all the HH residuals...
csvwrite('C:\G-LUM Finale\Calibration\Output\ResHH.txt',eps);
'Finished Calibrating Residential Location Choice'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EMPLOC%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% The path of input& output file must be specified ..
% The path of input& output file must be specified ..

clear all;

HH_lag = load('C:\G-LUM Finale\Calibration\Input\HH_lag.txt'); % Inputfile path must be specified
EMPactual = load('C:\G-LUM Finale\Calibration\Input\EMP_base.txt');% Inputfile path must be specified
EMP_lag = load('C:\G-LUM Finale\Calibration\Input\EMP_lag.txt');% Inputfile path must be specified

%Dimensions of the HH & EMP matrices
dimZone=size(HH_lag,1); 
dimHH=size(HH_lag,2); 
dimEMP=size(EMP_lag,2);
empclass=dimEMP;

%Fill in 0.0001 when HH or EMP is 0
for i=1:dimZone
    for j=1:dimHH
        if HH_lag(i,j) < 0.0001
            HH_lag(i,j)=0.0001;
        end;
    end;
    
    for k=1:dimEMP
        if EMP_lag(i,k) < 0.0001
            EMP_lag(i,k)=0.0001;
        end;
        if EMPactual(i,k) < 0.0001
            EMPactual(i,k)=0.0001;
        end;
    end;
end;

if empclass>0 EMP_lag_1=EMP_lag(:,1); end
if empclass>1 EMP_lag_2=EMP_lag(:,2); end
if empclass>2 EMP_lag_3=EMP_lag(:,3); end
if empclass>3 EMP_lag_4=EMP_lag(:,4); end
if empclass>4 EMP_lag_5=EMP_lag(:,5); end
if empclass>5 EMP_lag_6=EMP_lag(:,6); end


if empclass>0 EMPactual_1=EMPactual(:,1); end
if empclass>1 EMPactual_2=EMPactual(:,2); end
if empclass>2 EMPactual_3=EMPactual(:,3); end
if empclass>3 EMPactual_4=EMPactual(:,4); end
if empclass>4 EMPactual_5=EMPactual(:,5); end
if empclass>5 EMPactual_6=EMPactual(:,6); end

%Calculate zonal total HH
sumHHzone=sum(HH_lag,2); %1*** by 1

%Zone size
ZoneSize = load('C:\G-LUM Finale\Calibration\Input\ZoneSize.txt');% Inputfile path must be specified

%Impedance 
c = load('C:\G-LUM Finale\Calibration\Input\TT.txt');% Inputfile path must be specified

options = optimset('MaxIter',1000000, 'MaxFunEvals', 1000000,'Tolx', 0.001,'TolFun',0.001);

%Initialization
%lambda=0.05; omega= 2.00 ; rho=-0.001;
%aa=0.5; bb= 0.05; 
%y0=vertcat(lambda,omega,rho,aa,bb);

if empclass>0
%z = @(y) NLLS_EMPLOC(y,EMP_lag_1,EMPactual_1,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
%[y,fval,exitflag]=fminsearch(z,y0,options);
%csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_1_fval.dat',fval); % Output file path must be specified
%csvwrite('C:\G-LUM Finale\Calibration\Output\EMPALLOC_1_parameters.dat',y); % output file path must be specified

%Initialization
%lambda=0.05; omega= 2.00 ; rho=-0.001;    %lambda=iv(5,1);omega=iv(5,2);rho=iv(5,3);aa=iv(5,4);bb=iv(5,5);
%aa=0.5; bb= 0.05   ; 
iv=load('C:\G-LUM Finale\Calibration\Input\initial_values.txt');
lambda=iv(7,1); omega=iv(7,2); rho=iv(7,3);
aa=iv(7,4); bb=iv(7,5); 
y0=vertcat(lambda,omega,rho,aa,bb);
z = @(y) NLLS_EMPLOC(y,EMP_lag_1,EMPactual_1,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
[y,fval,exitflag]=fminsearch(z,y0,options);

csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_1_fval.dat',fval); % Output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_1_parameters.dat',y); % output file path must be specified
end

if empclass>1
%Initialization
%lambda=0.05; omega= 2.00 ; rho=-0.001;
%aa=0.5; bb= 0.05; 
lambda=iv(8,1); omega=iv(8,2); rho=iv(8,3);
aa=iv(8,4); bb=iv(8,5); 
y0=vertcat(lambda,omega,rho,aa,bb);

z = @(y)NLLS_EMPLOC(y,EMP_lag_2,EMPactual_2,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
[y,fval,exitflag]=fminsearch(z,y0,options);
end

if empclass>2
csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_2_fval.dat',fval);% Output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_2_parameters.dat',y); % output file path must be specified

%Initialization
%lambda=0.05; omega= 2.00 ; rho=-0.001;
%aa=0.5; bb= 0.05; 
lambda=iv(9,1); omega=iv(9,2); rho=iv(9,3);
aa=iv(9,4); bb=iv(9,5); 
y0=vertcat(lambda,omega,rho,aa,bb);

z = @(y)NLLS_EMPLOC(y,EMP_lag_3,EMPactual_3,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
[y,fval,exitflag]=fminsearch(z,y0,options);

csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_3_fval.dat',fval);% Output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_3_parameters.dat',y); % Output file path must be specified
end
if empclass>3
%Initialization
%lambda=0.05; omega= 2.00 ; rho=-0.001;
%aa=0.5; bb= 0.05; 
lambda=iv(10,1); omega=iv(10,2); rho=iv(10,3);
aa=iv(10,4); bb=iv(10,5); 

y0=vertcat(lambda,omega,rho,aa,bb);

z = @(y)NLLS_EMPLOC(y,EMP_lag_4,EMPactual_4,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
[y,fval,exitflag]=fminsearch(z,y0,options);

csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_4_fval.dat',fval);% Output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_4_parameters.dat',y); % Output file path must be specified
end
if empclass>4
% Initialization
% lambda=0.05; omega= 2.00 ; rho=-0.001;% aa=0.5; bb= 0.05; 
lambda=iv(11,1); omega=iv(11,2); rho=iv(11,3); aa=iv(11,4); bb=iv(11,5); 

y0=vertcat(lambda,omega,rho,aa,bb);
z = @(y)NLLS_EMPLOC(y,EMP_lag_5,EMPactual_5,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
[y,fval,exitflag]=fminsearch(z,y0,options);

csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_5_fval.dat',fval);% Output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_5_parameters.dat',y); % output file path must be specified
end

if empclass>5
lambda=iv(12,1); omega=iv(12,2); rho=iv(12,3); aa=iv(12,4); bb=iv(12,5); 
y0=vertcat(lambda,omega,rho,aa,bb);

z = @(y)NLLS_EMPLOC(y,EMP_lag_6,EMPactual_6,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
[y,fval,exitflag]=fminsearch(z,y0,options);

csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_6_fval.dat',fval);% Output file path must be specified
csvwrite('C:\G-LUM Finale\Calibration\Output\EMPLOC_6_parameters.dat',y); % output file path must be specified
end
% Calculation of residuals 

if empclass>0
    % Employment Category 1 residual computation
    y0 = load('C:\G-LUM Finale\Calibration\Output\EMPLOC_1_parameters.dat');
    z= NLLS_EMPLOC_SIM(y0,EMP_lag_1,EMPactual_1,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
    eps_Emp1 = load('C:\G-LUM Finale\Calibration\Output\epsilon_emp.dat');
end
if empclass>1
    % Employment Category 2 residual computation
    y0 = load('C:\G-LUM Finale\Calibration\Output\EMPLOC_2_parameters.dat');
    z= NLLS_EMPLOC_SIM(y0,EMP_lag_2,EMPactual_2,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
    eps_Emp2 = load('C:\G-LUM Finale\Calibration\Output\epsilon_emp.dat');
end
if empclass>2
    % Employment Category 3 residual computation
    y0 = load('C:\G-LUM Finale\Calibration\Output\EMPLOC_3_parameters.dat');
    z= NLLS_EMPLOC_SIM(y0,EMP_lag_3,EMPactual_3,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
    eps_Emp3 = load('C:\G-LUM Finale\Calibration\Output\epsilon_emp.dat');
end
if empclass>3
    % Employment Category 4 residual computation
    y0 = load('C:\G-LUM Finale\Calibration\Output\EMPLOC_4_parameters.dat');
    z= NLLS_EMPLOC_SIM(y0,EMP_lag_4,EMPactual_4,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
    eps_Emp4 = load('C:\G-LUM Finale\Calibration\Output\epsilon_emp.dat');
end
if empclass>4
    % Employment Category 5 residual computation
    y0 = load('C:\G-LUM Finale\Calibration\Output\EMPLOC_5_parameters.dat');
    z= NLLS_EMPLOC_SIM(y0,EMP_lag_5,EMPactual_5,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
    eps_Emp5 = load('C:\G-LUM Finale\Calibration\Output\epsilon_emp.dat');
end
if empclass>5
    % Employment Category 6 residual computation
    y0 = load('C:\G-LUM Finale\Calibration\Output\EMPLOC_6_parameters.dat');
    z= NLLS_EMPLOC_SIM(y0,EMP_lag_6,EMPactual_6,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);
    eps_Emp6 = load('C:\G-LUM Finale\Calibration\Output\epsilon_emp.dat');
end

% Integrating all the residuals ..
%eps_Emp = [eps_Emp1';eps_Emp2';eps_Emp3';eps_Emp4';eps_Emp5';eps_Emp6']';
%%Siva Modified
if empclass>0 eps_Emp = [eps_Emp1';]'; end 
if empclass>1 eps_Emp = [eps_Emp1';eps_Emp2']'; end 
if empclass>2 eps_Emp = [eps_Emp1';eps_Emp2';eps_Emp3']'; end 
if empclass>3 eps_Emp = [eps_Emp1';eps_Emp2';eps_Emp3';eps_Emp4']'; end 
if empclass>4 eps_Emp = [eps_Emp1';eps_Emp2';eps_Emp3';eps_Emp4';eps_Emp5']'; end 
if empclass>5 eps_Emp = [eps_Emp1';eps_Emp2';eps_Emp3';eps_Emp4';eps_Emp5';eps_Emp6']'; end 
csvwrite('C:\G-LUM Finale\Calibration\Output\ResEMP.txt',eps_Emp);
'Finished Calibrating Employment Location Choice'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%CALIBRATION LUDENSITY%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
clear all;
c = load('C:\G-LUM Finale\Prediction\Input\TT.txt'); % Travel time..% Input path must be specified
HHactual = load('C:\G-LUM Finale\Calibration\Input\HH_base.txt');%Input path must be specified
EMPactual = load('C:\G-LUM Finale\Calibration\Input\EMP_base.txt');%Input path must be specified


%Dimensions of the HH & EMP matrices
dimZone=size(HHactual,1); %It should be 1***.
dimHH=size(HHactual,2); %It should be 6.
dimEMP=size(EMPactual,2); %It should be 3.
hhclass=dimHH;
emclass=dimEMP;

for i=1:dimZone
    for j=1:dimHH
        if HHactual(i,j) < 0.0001
            HHactual(i,j)=0.0001;
        end;
    end;
    
    for k=1:dimEMP
        if EMPactual(i,k) < 0.0001
            EMPactual(i,k)=0.0001;
        end;
    end;
end;

%Calculate HH by type
if hhclass>0 HH_1=HHactual(:,1); end
if hhclass>1 HH_2=HHactual(:,2); end
if hhclass>2 HH_3=HHactual(:,3); end
if hhclass>3 HH_4=HHactual(:,4); end
if hhclass>4 HH_5=HHactual(:,5); end
if hhclass>5 HH_6=HHactual(:,6); end

%Calculate EMP by BAS & COMM
BASzone=EMPactual(:,1);
COMMzone=EMPactual(:,2)+EMPactual(:,3); 

%Calculate zonal total HH
sumHHzone=sum(HHactual,2); %1*** by 1
sumEMPzone=sum(EMPactual,2); %1*** by 1

%All necessary land in each zone
landrP = load('C:\G-LUM Finale\Calibration\Input\landr.txt'); %residential..Input path must be specified
landbP = load('C:\G-LUM Finale\Calibration\Input\landb.txt'); %basic..Input path must be specified
landeP = load('C:\G-LUM Finale\Calibration\Input\landu.txt'); %unusable land..Input path must be specified
landuP = load('C:\G-LUM Finale\Calibration\Input\landd.txt'); %undeveloped ..developable..Input path must be specified
landcP = load('C:\G-LUM Finale\Calibration\Input\landc.txt'); %commerical..Input path must be specified
landsP = load('C:\G-LUM Finale\Calibration\Input\lands.txt'); %Streets & Highways..Input path must be specified

for i=1:dimZone
    if(landrP(i) < 0.0001) landrP(i) = 0.0001; end;
    if(landbP(i) < 0.0001) landbP(i) = 0.0001; end;
    if(landeP(i) < 0.0001) landeP(i) = 0.0001; end;
    if(landuP(i) < 0.0001) landuP(i) = 0.0001; end;
    if(landcP(i) < 0.0001) landcP(i) = 0.0001; end;
    if(landsP(i) < 0.0001) landsP(i) = 0.0001; end;
end;

ZoneSize = load('C:\G-LUM Finale\Calibration\Input\ZoneSize.txt');%Input path must be specified

landdable = ZoneSize-landeP; % Developable land ..
landd = ZoneSize - landuP - landeP; % Amount of developed land..
landb = landbP; % Land for basic employment..
landc = landcP;% Land for commercial employment..
landr = landrP;

%Initialization

iv=load('C:\G-LUM Finale\Calibration\Input\initial_values.txt');
kk=iv(13,1).*ones(1,5+dimHH); gg=iv(13,2).*ones(1,5+dimHH); pp=iv(13,3)*ones(1,4+dimHH); 
k0=kk'; g0=gg'; p0=pp';
options = optimset('MaxIter',1000000, 'MaxFunEvals', 1000000,'Tolx', 0.001,'TolFun',0.001);

z1 = @(k) NLLS_LUDENSITY_K(k,landr,landdable,landd,landb,landc,sumHHzone,HHactual,dimZone,sumEMPzone,BASzone,COMMzone);
[k,fval,exitflag]=fminsearch(z1,k0,options);
csvwrite('C:\G-LUM Finale\Calibration\output\LUDENSITY_k_fval.dat',fval);%Output path must be specified
csvwrite('C:\G-LUM Finale\Calibration\output\LUDENSITY_k_param.dat',k);%Output path must be specified
%csvwrite('C:\G-LUM Finale\Prediction\Output\test1.dat',HHactual(7,2));
z2 = @(g) NLLS_LUDENSITY_G(g,landr,landdable,landd,landb,sumEMPzone,BASzone,COMMzone,dimZone,HHactual,sumHHzone);
[g,fval,exitflag]=fminsearch(z2,g0,options);
%csvwrite('C:\G-LUM Finale\Prediction\Output\test.dat',g);

csvwrite('C:\G-LUM Finale\Calibration\output\LUDENSITY_g_fval.dat',fval);%Output path must be specified
csvwrite('C:\G-LUM Finale\Calibration\output\LUDENSITY_g_param.dat',g);%Output path must be specified

z3 = @(p) NLLS_LUDENSITY_P(p,landr,landdable,landd,landc,sumEMPzone,BASzone,COMMzone,dimZone,HHactual,sumHHzone);
[p,fval,exitflag]=fminsearch(z3,p0,options);

csvwrite('C:\G-LUM Finale\Calibration\output\LUDENSITY_p_fval.dat',fval);%Output path must be specified
csvwrite('C:\G-LUM Finale\Calibration\output\LUDENSITY_p_param.dat',p);%Output path must be specified

% residual computations..
k0 = load('C:\G-LUM Finale\Calibration\output\LUDENSITY_k_param.dat');
z1 = NLLS_LUDENSITY_K_SIM(k0,landr,landdable,landd,landb,landc,sumHHzone,HHactual,dimZone,sumEMPzone,BASzone,COMMzone);
eps_K = load('C:\G-LUM Finale\Calibration\output\epsilon_k.dat');

g0 = load('C:\G-LUM Finale\Calibration\output\LUDENSITY_g_param.dat');
z2 = NLLS_LUDENSITY_G_SIM(g0,landr,landdable,landd,landb,sumEMPzone,BASzone,COMMzone,dimZone,HHactual,sumHHzone);
eps_G = load('C:\G-LUM Finale\Calibration\output\epsilon_g.dat');

p0 = load('C:\G-LUM Finale\Calibration\output\LUDENSITY_p_param.dat');
z3 = NLLS_LUDENSITY_P_SIM(p0,landr,landdable,landd,landc,sumEMPzone,BASzone,COMMzone,dimZone,HHactual,sumHHzone);
eps_P = load('C:\G-LUM Finale\Calibration\output\epsilon_p.dat');

% Intergrating all the residuals ..
eps_LU = [eps_K';eps_G';eps_P']';
csvwrite('C:\G-LUM Finale\Calibration\output\ResLU.txt',eps_LU);
'Finished Calibrating LUDENSITY model'

%*********************************************PREDICTION****************************************
% --- Executes on button press in pushbutton4.
function pushbutton4_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton4 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% The path for input & output file must be speficied
% The path for input & output file must be speficied
% The path for input & output file must be speficied
clear all;
for ii=1:5:5
        UncertainReal=ii;
        jj=0;
        jj=jj+1;
        %Input
        HHbase_original = load('C:\G-LUM Finale\Prediction\Input\HH_base.txt'); % Input path must be specified
        EMPbase_original = load('C:\G-LUM Finale\Prediction\Input\EMP_base.txt');% Input path must be specified

        %Dimensions of the HH & EMP matrices
        dimZone=size(HHbase_original,1); 
        dimHH=size(HHbase_original,2);
        dimEMP=size(EMPbase_original,2);
        hhclass=dimHH;
        empclass=dimEMP;
        % Unemployment rate is assumed to be zero
        unemp=zeros(1,dimEMP);

        %Base year land use data
        landrP = load('C:\G-LUM Finale\Prediction\Input\landr.txt'); %residential ..% Input path must be specified
        landbP = load('C:\G-LUM Finale\Prediction\Input\landb.txt'); %basic..% Input path must be specified
        landeP = load('C:\G-LUM Finale\Prediction\Input\landu.txt'); %unusable land..% Input path must be specified
        landuP = load('C:\G-LUM Finale\Prediction\Input\landd.txt'); %undeveloped or developable..% Input path must be specified
        landcP = load('C:\G-LUM Finale\Prediction\Input\landc.txt'); %commerical..% Input path must be specified
        landsP = load('C:\G-LUM Finale\Prediction\Input\lands.txt'); %Streets & Highways..% Input path must be specified

        for i=1:dimZone
            if(landrP(i) < 0.001) landrP(i) = 0.001;end;
            if(landbP(i) < 0.001) landbP(i) = 0.001;end;
            if(landeP(i) < 0.001) landeP(i) = 0.001;end;
            if(landuP(i) < 0.001) landuP(i) = 0.001;end;
            if(landcP(i) < 0.001) landcP(i) = 0.001;end;
            if(landsP(i) < 0.001) landsP(i) = 0.001;end;
        end;

        landv = landuP; 
        landr = landrP;
        landdRatio = (landrP + landbP  + landcP + landsP )./( landrP + landbP  + landuP + landcP+ landsP);  % Proportion of developable land developed
        ZoneSize = load('C:\G-LUM Finale\Prediction\Input\ZoneSize.txt');% Input path must be specified

        landdable = landrP + landbP  + landuP + landcP+ landsP; % Developable land
        landd = landrP + landbP + landcP+ landsP; % Amount of developed land
        landb = landbP; % Land for basic employment
        landc = landcP;% Land for commercial employment

        % Impedance data
        c = load('C:\G-LUM Finale\Prediction\Input\TT.txt'); % Travel time..% Input path must be specified
        % RESLOC parameters
        % Parameters obtained by Calibration .....

        if hhclass>0
            HH_Parameter_temp1 = load('C:\G-LUM Finale\Prediction\Input\RESLOC_1_parameters.dat');
            HH_Parameter = [HH_Parameter_temp1'];
        end
        if hhclass>1 
            HH_Parameter_temp2 = load('C:\G-LUM Finale\Prediction\Input\RESLOC_2_parameters.dat');
            HH_Parameter = [HH_Parameter_temp1';HH_Parameter_temp2'];
        end
        if hhclass>2
            HH_Parameter_temp3 = load('C:\G-LUM Finale\Prediction\Input\RESLOC_3_parameters.dat');
            HH_Parameter = [HH_Parameter_temp1';HH_Parameter_temp2';HH_Parameter_temp3']    ;
        end
        if hhclass>3
            HH_Parameter_temp4 = load('C:\G-LUM Finale\Prediction\Input\RESLOC_4_parameters.dat');
            HH_Parameter = [HH_Parameter_temp1';HH_Parameter_temp2';HH_Parameter_temp3';HH_Parameter_temp4'];
        end
        if hhclass>4
            HH_Parameter_temp5 = load('C:\G-LUM Finale\Prediction\Input\RESLOC_5_parameters.dat');
            HH_Parameter = [HH_Parameter_temp1';HH_Parameter_temp2';HH_Parameter_temp3';HH_Parameter_temp4';HH_Parameter_temp5'];
        end
        if hhclass>5
            HH_Parameter_temp6 = load('C:\G-LUM Finale\Prediction\Input\RESLOC_6_parameters.dat');
            HH_Parameter = [HH_Parameter_temp1';HH_Parameter_temp2';HH_Parameter_temp3';HH_Parameter_temp4';HH_Parameter_temp5';HH_Parameter_temp6'];
        end
             
        % EMPLOC paramters
        % Parameter obtained by Calibration 

        if empclass>0 EMP_Parameter_temp1 = load('C:\G-LUM Finale\Prediction\Input\EMPLOC_1_parameters.dat'); end
        if empclass>1 EMP_Parameter_temp2 = load('C:\G-LUM Finale\Prediction\Input\EMPLOC_2_parameters.dat'); end 
        if empclass>2 EMP_Parameter_temp3 = load('C:\G-LUM Finale\Prediction\Input\EMPLOC_3_parameters.dat'); end 
        if empclass>3 EMP_Parameter_temp4 = load('C:\G-LUM Finale\Prediction\Input\EMPLOC_4_parameters.dat'); end
        if empclass>4 EMP_Parameter_temp5 = load('C:\G-LUM Finale\Prediction\Input\EMPLOC_5_parameters.dat'); end
        if empclass>5 EMP_Parameter_temp6 = load('C:\G-LUM Finale\Prediction\Input\EMPLOC_6_parameters.dat'); end
        
        if empclass>0 EMP_Parameter = [EMP_Parameter_temp1']; end
        if empclass>1 EMP_Parameter = [EMP_Parameter_temp1';EMP_Parameter_temp2']; end 
        if empclass>2 EMP_Parameter = [EMP_Parameter_temp1';EMP_Parameter_temp2';EMP_Parameter_temp3']; end 
        if empclass>3 EMP_Parameter = [EMP_Parameter_temp1';EMP_Parameter_temp2';EMP_Parameter_temp3';EMP_Parameter_temp4']; end
        if empclass>4 EMP_Parameter = [EMP_Parameter_temp1';EMP_Parameter_temp2';EMP_Parameter_temp3';EMP_Parameter_temp4';EMP_Parameter_temp5']; end
        if empclass>5 EMP_Parameter = [EMP_Parameter_temp1';EMP_Parameter_temp2';EMP_Parameter_temp3';EMP_Parameter_temp4';EMP_Parameter_temp5';EMP_Parameter_temp6']; end


        % Loading the residuals ...
        eps_HH = load('C:\G-LUM Finale\Prediction\Input\ResHH.txt'); % residuals of the households
        eps_EMP = load('C:\G-LUM Finale\Prediction\Input\ResEMP.txt'); % residuals of the employment
        eps_LU = load('C:\G-LUM Finale\Prediction\Input\ResLU.txt'); % residuals of the land cover



        eta=HH_Parameter(:,1)';
        alpha=HH_Parameter(:,2)';
        beta=HH_Parameter(:,3)';
        qq=HH_Parameter(:,4)';
        rr=HH_Parameter(:,5)';
        ss=HH_Parameter(:,6)';
        bnn=HH_Parameter(:,7:6+dimHH)';

        lambda=EMP_Parameter(:,1)';
        omega=EMP_Parameter(:,2)';
        rho=EMP_Parameter(:,3)';
        aa=EMP_Parameter(:,4)';
        bb=EMP_Parameter(:,5)';

        % LUDENSITY parameters ....
        % Paramters obtained by Calibration..
        k = load('C:\G-LUM Finale\Prediction\Input\LUDENSITY_k_param.dat');
        % Paramters obtained by Calibration..
        g = load('C:\G-LUM Finale\Prediction\Input\LUDENSITY_g_param.dat');
        % Paramters obtained by Calibration..
        p = load('C:\G-LUM Finale\Prediction\Input\LUDENSITY_p_param.dat');

        %Initialization
    %   k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);  k7=k(8,1); k8=k(9,1); k9=k(10,1); k10=k(11,1); 
     %  g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);g6=g(7,1);g7=g(8,1);g8=g(9,1);g9=g(10,1);
    %   p0=p(1,1); p1=p(2,1); p2=p(3,1); p3=p(4,1); p4=p(5,1); p5=p(6,1); p6=p(7,1); p7=p(8,1); p8=p(9,1); p9=p(10,1);

if dimHH==1
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1);
   g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); 
   p0=p(1,1); p1=p(2,1); p2=p(3,1); p3=p(4,1); p4=p(5,1);
elseif dimHH==2
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);
    p0=p(1,1); p1=p(2,1); p2=p(3,1); p3=p(4,1); p4=p(5,1); p5=p(6,1);
elseif dimHH==3
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);  k7=k(8,1);
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);g6=g(7,1);
    p0=p(1,1); p1=p(2,1); p2=p(3,1); p3=p(4,1); p4=p(5,1); p5=p(6,1); p6=p(7,1);
elseif dimHH==4
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);  k7=k(8,1); k8=k(9,1); 
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);g6=g(7,1);g7=g(8,1);
    p0=p(1,1); p1=p(2,1); p2=p(3,1); p3=p(4,1); p4=p(5,1); p5=p(6,1); p6=p(7,1); p7=p(8,1); 
elseif dimHH==5
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);  k7=k(8,1); k8=k(9,1); k9=k(10,1);
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);g6=g(7,1);g7=g(8,1);g8=g(9,1);
    p0=p(1,1); p1=p(2,1); p2=p(3,1); p3=p(4,1); p4=p(5,1); p5=p(6,1); p6=p(7,1); p7=p(8,1); p8=p(9,1);
elseif dimHH==6
    k0=k(1,1); k1=k(2,1); k2=k(3,1); k3=k(4,1); k4=k(5,1); k5=k(6,1); k6=k(7,1);  k7=k(8,1); k8=k(9,1); k9=k(10,1); k10=k(11,1); 
    g0=g(1,1); g1=g(2,1); g2=g(3,1); g3=g(4,1); g4=g(5,1); g5=g(6,1);g6=g(7,1);g7=g(8,1);g8=g(9,1);g9=g(10,1);
    p0=p(1,1); p1=p(2,1); p2=p(3,1); p3=p(4,1); p4=p(5,1); p5=p(6,1); p6=p(7,1); p7=p(8,1); p8=p(9,1); p9=p(10,1);
end 
       
       
       
       
        % Regional control totals
        % Control Totals for 2010,2015,2020,2025,2030....
        EMPprojectCT = load('C:\G-LUM Finale\Prediction\Input\EmpProj.txt');
        HHprojectCT = load('C:\G-LUM Finale\Prediction\Input\HHProj.txt');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FUTURE PREDICTIONS ....................................

        % Multi-period predictions ....
           dimTimePd=load('C:\G-LUM Finale\Prediction\Output\years.txt');
           for count = 1:dimTimePd
                %    ProjectionYear = 0+ (5*count)  %(Jen) changed from 2005 + (5*count), to generalize
                % Seeding 
                seed=1.0.*ones(dimZone,1);

                HHbase=HHbase_original; 
                EMPbase=EMPbase_original;
                HHbase_seed=zeros(dimZone,dimHH); 
                EMPbase_seed=zeros(dimZone,dimEMP);

                for i=1:dimZone;
                    for n=1:dimHH;
                        if HHbase_original(i,n) < seed(i,1);
                           HHbase(i,n)=HHbase(i,n)+seed(i,1);
                           HHbase_seed(i,n)=1;
                        end;
                    end;

                    for k=1:dimEMP;
                        if EMPbase_original(i,k) < seed(i,1);
                           EMPbase(i,k)=EMPbase(i,k)+seed(i,1);
                           EMPbase_seed(i,k)=1;
                        end;
                    end;    

                end;

        %Calculate zonal total HH
        sumHHzone=sum(HHbase,2); %1*** by 1


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%EMPLOC%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Calculation in EMPLOC 
%Calculate A
        Ainv=zeros(dimZone,dimEMP);
                for i=1:dimZone;
                    for k=1:dimEMP;
                        for j=1:dimZone;
                            Ainv(i,k)=Ainv(i,k)+EMPbase(j,k).^aa(1,k).*ZoneSize(j,1).^bb(1,k).*c(i,j).^omega(1,k).*exp(rho(1,k)*c(i,j));
                        end;
                        A(i,k)=1./Ainv(i,k);
                    end;
                end;

                %Calculate W_EMPLOC
                for j=1:dimZone;
                    for k=1:dimEMP;
                        W_EMPLOC(j,k)=EMPbase(j,k).^aa(1,k).*ZoneSize(j,1).^bb(1,k);     
                    end;
                end;
%Calculate modeled EMP
                EMPsum=zeros(dimZone,dimEMP);
                for j=1:dimZone;
                    for k=1:dimEMP;
                        for i=1:dimZone;
                            EMPsum(j,k)=EMPsum(j,k)+ (sumHHzone(i,1).*A(i,k).*W_EMPLOC(j,k).*c(i,j).^omega(1,k)*exp(rho(1,k).*c(i,j)));
                        end;
                        EMPproject_withseed(j,k)=(lambda(1,k).*EMPsum(j,k))+((1-lambda(1,k))*EMPbase(j,k));
                    end;
                end;

                %Deduct the seeds
                EMPproject_withoutseed = EMPproject_withseed;
                for i=1:dimZone;
                    for k=1:dimEMP;
                        if EMPbase_seed(i,k) == 1;
                           EMPproject_withoutseed(i,k)=max(0,EMPproject_withoutseed(i,k)-seed(i,1));
                        end;
                    end;        
                end;


                %Normalization WITH PROJECTED VALUES ......
                EMPproject_withoutseed_total=sum(EMPproject_withoutseed,1); %1 by 3

                 for i=1:dimZone;
                     for k=1:dimEMP;
                         EMPproject(i,k)=EMPproject_withoutseed(i,k).*(EMPprojectCT(count,k)./EMPproject_withoutseed_total(1,k));
                     end;
                 end;
        
 % Adding the residuals .....
                 if(count == 1) temp_EMP = EMPproject + eps_EMP; end;
                 if(count == 2) temp_EMP = EMPproject + (0.75*eps_EMP);end;
                 if(count == 3) temp_EMP = EMPproject + (0.5*eps_EMP);end;
                 if(count == 4) temp_EMP = EMPproject + (0.25*eps_EMP); end;
                 if(count == 5) temp_EMP = EMPproject ; end;

                  for i=1:dimZone;
                     for k=1:dimEMP;
                         if(temp_EMP(i,k) > 0) EMPproject(i,k) =  temp_EMP(i,k);end;
                         if(EMPproject(i,k) < 0) EMPproject(i,k) = 0; end;
                     end;
                  end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%END OF EMPLOC%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%RESLOC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Calculation in RESLOC 
%Calculate "a" (ratio of HH to EMP)

                sumHHtype=sum(HHbase_original,1); %1 by dimHH
                sumEMPtype=sum(EMPbase_original,1); %1 by dimHH
                a=(ones(1,dimEMP)./sumEMPtype)'*sumHHtype; %4 by 4

                %Calculate Q
                Q=zeros(dimZone,dimHH);
                for j=1:dimZone;
                    for n=1:dimHH;
                        for k=1:dimEMP;
                            Q(j,n)=Q(j,n)+((a(k,n).*EMPproject(j,k)./(1-unemp(1,k))));
                        end;
                    end;
                end;

                %Calculate W_RESLOC
                pipi=ones(dimZone,dimHH); %Initialize the multi-time part
                for i=1:dimZone;
                    for n=1:dimHH;
                        for nn=1:dimHH;
                            pipi(i,n)=pipi(i,n).*((1+HHbase(i,nn)./sumHHzone(i,1)).^bnn(nn,n));
                        end;

                        W_RESLOC(i,n)=landv(i,1).^qq(1,n)*(1+landdRatio(i,1)).^rr(1,n)*landr(i,1).^ss(1,n)*pipi(i,n);     
                    end;
                end;

                %Calculate B
                Binv=zeros(dimZone,dimHH);
                for j=1:dimZone;
                    for n=1:dimHH;
                        for i=1:dimZone;
                            Binv(j,n)=Binv(j,n)+(W_RESLOC(i,n).*c(i,j).^alpha(1,n).*exp(beta(1,n)*c(i,j)));
                        end;
                        B(j,n)=1./Binv(j,n);
                    end;
                end;

        %Calculate modeled HH
        HHsum=zeros(dimZone,dimHH);
                for i=1:dimZone;
                    for n=1:dimHH;
                        for j=1:dimZone;
                            value(j)=B(j,n);
                            value1(j)=W_RESLOC(i,n).*c(i,j).^alpha(1,n).*exp(beta(1,n)*c(i,j));
                             HHsum(i,n)=HHsum(i,n)+ (Q(j,n).*value(j)).*value1(j);
                            %HHsum(i,n)=HHsum(i,n)+ (Q(j,n).*B(j,n).*W_RESLOC(i,n).*c(i,j).^alpha(1,n).*exp(beta(1,n)*c(i,j)));
                           if i==31
                               if n==3
                            HHcheck(j)=HHsum(i,n);
                           % value(j)= (Q(j,n).*B(j,n).*W_RESLOC(i,n).*c(i,j).^alpha(1,n).*exp(beta(1,n)*c(i,j)));
                           % value1(j)=W_RESLOC(i,n).*c(i,j).^alpha(1,n).*exp(beta(1,n)*c(i,j));
                                end;
                            end;
                        end;
                        HHproject_withseed(i,n)=(eta(1,n).*HHsum(i,n))+((1-eta(1,n)).*HHbase(i,n));
                    end;
                end;

                %Deduct the seeds
                HHproject_withoutseed=HHproject_withseed;
                for i=1:dimZone;
                    for n=1:dimHH;
                        if HHbase_seed(i,n) == 1;
                           HHproject_withoutseed(i,n)=max(0,HHproject_withoutseed(i,n)-seed(i,1));
                        end;
                    end;        
                end;

                %Normalization
                 HHproject_withoutseed_total=sum(HHproject_withoutseed,1); %1 by 6

                 for i=1:dimZone;
                     for n=1:dimHH;
                         HHproject(i,n)= HHproject_withoutseed(i,n).*(HHprojectCT(count,n)./HHproject_withoutseed_total(1,n));
                     end;
                 end;
csvwrite('C:\G-LUM Finale\Prediction\Output\check.dat',HHcheck'); 
%csvwrite('C:\G-LUM Finale\Prediction\Output\check1.dat',value');
%csvwrite('C:\G-LUM Finale\Prediction\Output\check2.dat',value1');
        str=input('count')
% Adding the residuals .....
                if(count == 1) temp_HH = HHproject + eps_HH; end;
                if(count == 2) temp_HH = HHproject + (0.75*eps_HH);end;
                if(count == 3) temp_HH = HHproject + (0.5*eps_HH);end;
                if(count == 4) temp_HH = HHproject + (0.25*eps_HH);end;
                if(count == 5) temp_HH = HHproject;end;


                for i=1:dimZone;
                     for n=1:dimHH;
                        if(temp_HH(i,n) > 0 )  HHproject(i,n) = temp_HH(i,n);end;
                        if(HHproject(i,n) < 0 ) HHproject(i,n) = 0; end;
                     end;
                end;


%%%%%%%%%%%%%%%%% END OF RESLOC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%% LUDENSITY %%%%%%%%%%%%%%%%%%%%%%%

%Calculation in LUDENSITY ....

                HHproject_LU = HHproject;
                EMPproject_LU = EMPproject;
                for i=1:dimZone
                    for n=1:dimHH
                        if HHproject_LU(i,n) < 0.0001
                           HHproject_LU(i,n)=0.0001;
                        end;
                    end;

                    for k=1:dimEMP
                        if EMPproject_LU(i,k) < 0.0001
                           HHproject_LU(i,k)=0.0001;
                        end;
                    end;   
                end;

                % HH & EMP in 2010
                HH_zone1 = HHproject_LU(1,:); 

                if hhclass>0        HH_1=HHproject_LU(:,1);         end
                if hhclass>1        HH_2=HHproject_LU(:,2);         end
                if hhclass>2        HH_3=HHproject_LU(:,3);         end
                if hhclass>3        HH_4=HHproject_LU(:,4);         end
                if hhclass>4        HH_5=HHproject_LU(:,5);         end
                if hhclass>5        HH_6=HHproject_LU(:,6);         end

                sumHHzoneproject=sum(HHproject_LU,2); %1*** by 1
                sumEMPzoneproject=sum(EMPproject_LU,2); %1*** by 1
                BASzoneproject=EMPproject_LU(:,1);
                COMMzoneproject=EMPproject_LU(:,2)+EMPproject_LU(:,3);

                for i=1:dimZone
                            if(landd(i) < 0.001) landd(i) = 0.001;end;
                            if(landdable(i) < 0.001) landdable(i) = 0.001;end;
                            if(sumHHzoneproject(i) < 0.001) sumHHzoneproject(i) = 0.001;end;
                            if(sumEMPzoneproject(i) < 0.001) sumEMPzoneproject(i) = 0.001;end;
                end;
                %Calculate Lr

                
                if hhclass==1
                Lr_unrestricted =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzoneproject./sumEMPzoneproject).^k3).*((COMMzoneproject./sumEMPzoneproject).^k4).*((HH_1./sumHHzoneproject).^k5));
                elseif hhclass==2
                Lr_unrestricted =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzoneproject./sumEMPzoneproject).^k3).*((COMMzoneproject./sumEMPzoneproject).^k4).*((HH_1./sumHHzoneproject).^k5).*((HH_2./sumHHzoneproject).^k6).*sumHHzoneproject) ;
                elseif hhclass==3
                Lr_unrestricted =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzoneproject./sumEMPzoneproject).^k3).*((COMMzoneproject./sumEMPzoneproject).^k4).*((HH_1./sumHHzoneproject).^k5).*((HH_2./sumHHzoneproject).^k6).*((HH_3./sumHHzoneproject).^k7).*sumHHzoneproject) ;
                elseif hhclass==4
                Lr_unrestricted =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzoneproject./sumEMPzoneproject).^k3).*((COMMzoneproject./sumEMPzoneproject).^k4).*((HH_1./sumHHzoneproject).^k5).*((HH_2./sumHHzoneproject).^k6).*((HH_3./sumHHzoneproject).^k7).*((HH_4./sumHHzoneproject).^k8).*sumHHzoneproject) ;
                elseif hhclass==5 
                Lr_unrestricted =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzoneproject./sumEMPzoneproject).^k3).*((COMMzoneproject./sumEMPzoneproject).^k4).*((HH_1./sumHHzoneproject).^k5).*((HH_2./sumHHzoneproject).^k6).*((HH_3./sumHHzoneproject).^k7).*((HH_4./sumHHzoneproject).^k8).*((HH_5./sumHHzoneproject).^k9).*sumHHzoneproject) ;
                elseif hhclass==6 
                Lr_unrestricted =(k0.*(landdable.^k1).*((landd./landdable).^k2).*((BASzoneproject./sumEMPzoneproject).^k3).*((COMMzoneproject./sumEMPzoneproject).^k4).*((HH_1./sumHHzoneproject).^k5).*((HH_2./sumHHzoneproject).^k6).*((HH_3./sumHHzoneproject).^k7).*((HH_4./sumHHzoneproject).^k8).*((HH_5./sumHHzoneproject).^k9).*((HH_6./sumHHzoneproject).^k10).*sumHHzoneproject) ;
                end  
                
                %Calculate Lb
                if hhclass==1
                Lb_unrestricted = g0.*(landd./landdable).^g1.*(BASzoneproject./sumEMPzoneproject).^g2.*((COMMzoneproject./sumEMPzoneproject).^g3).*((HH_1./sumHHzoneproject).^g4).*BASzoneproject;
                elseif hhclass==2
                Lb_unrestricted = g0.*(landd./landdable).^g1.*(BASzoneproject./sumEMPzoneproject).^g2.*((COMMzoneproject./sumEMPzoneproject).^g3).*((HH_1./sumHHzoneproject).^g4).*((HH_2./sumHHzoneproject).^g5).*BASzoneproject;
                elseif hhclass==3
                Lb_unrestricted = g0.*(landd./landdable).^g1.*(BASzoneproject./sumEMPzoneproject).^g2.*((COMMzoneproject./sumEMPzoneproject).^g3).*((HH_1./sumHHzoneproject).^g4).*((HH_2./sumHHzoneproject).^g5).*((HH_3./sumHHzoneproject).^g6).*BASzoneproject;
                elseif hhclass==4
                Lb_unrestricted = g0.*(landd./landdable).^g1.*(BASzoneproject./sumEMPzoneproject).^g2.*((COMMzoneproject./sumEMPzoneproject).^g3).*((HH_1./sumHHzoneproject).^g4).*((HH_2./sumHHzoneproject).^g5).*((HH_3./sumHHzoneproject).^g6).*((HH_4./sumHHzoneproject).^g7).*BASzoneproject;
                elseif hhclass==5
                Lb_unrestricted = g0.*(landd./landdable).^g1.*(BASzoneproject./sumEMPzoneproject).^g2.*((COMMzoneproject./sumEMPzoneproject).^g3).*((HH_1./sumHHzoneproject).^g4).*((HH_2./sumHHzoneproject).^g5).*((HH_3./sumHHzoneproject).^g6).*((HH_4./sumHHzoneproject).^g7).*((HH_5./sumHHzoneproject).^g8).*BASzoneproject;
                elseif hhclass==6
                Lb_unrestricted = g0.*(landd./landdable).^g1.*(BASzoneproject./sumEMPzoneproject).^g2.*((COMMzoneproject./sumEMPzoneproject).^g3).*((HH_1./sumHHzoneproject).^g4).*((HH_2./sumHHzoneproject).^g5).*((HH_3./sumHHzoneproject).^g6).*((HH_4./sumHHzoneproject).^g7).*((HH_5./sumHHzoneproject).^g8).*((HH_6./sumHHzoneproject).^g9).*BASzoneproject;
                end 
                %Calculate Lc
                
                if hhclass==1
                Lc_unresricted = p0.*(landd./landdable).^p1.*(BASzoneproject./sumEMPzoneproject).^p2.*((COMMzoneproject./sumEMPzoneproject).^p3).*((HH_1./sumHHzoneproject).^p4).*COMMzoneproject;
                elseif hhclass==2
                Lc_unrestricted = p0.*(landd./landdable).^p1.*(BASzoneproject./sumEMPzoneproject).^p2.*((COMMzoneproject./sumEMPzoneproject).^p3).*((HH_1./sumHHzoneproject).^p4).*((HH_2./sumHHzoneproject).^p5).*COMMzoneproject;
                elseif hhclass==3
                Lc_unrestricted = p0.*(landd./landdable).^p1.*(BASzoneproject./sumEMPzoneproject).^p2.*((COMMzoneproject./sumEMPzoneproject).^p3).*((HH_1./sumHHzoneproject).^p4).*((HH_2./sumHHzoneproject).^p5).*((HH_3./sumHHzoneproject).^p6).*COMMzoneproject;
                elseif hhclass==4
                Lc_unrestricted = p0.*(landd./landdable).^p1.*(BASzoneproject./sumEMPzoneproject).^p2.*((COMMzoneproject./sumEMPzoneproject).^p3).*((HH_1./sumHHzoneproject).^p4).*((HH_2./sumHHzoneproject).^p5).*((HH_3./sumHHzoneproject).^p6).*((HH_4./sumHHzoneproject).^p7).*COMMzoneproject;
                elseif hhclass==5
                Lc_unrestricted = p0.*(landd./landdable).^p1.*(BASzoneproject./sumEMPzoneproject).^p2.*((COMMzoneproject./sumEMPzoneproject).^p3).*((HH_1./sumHHzoneproject).^p4).*((HH_2./sumHHzoneproject).^p5).*((HH_3./sumHHzoneproject).^p6).*((HH_4./sumHHzoneproject).^p7).*((HH_5./sumHHzoneproject).^p8).*COMMzoneproject;
                elseif hhclass==6
                Lc_unrestricted = p0.*(landd./landdable).^p1.*(BASzoneproject./sumEMPzoneproject).^p2.*((COMMzoneproject./sumEMPzoneproject).^p3).*((HH_1./sumHHzoneproject).^p4).*((HH_2./sumHHzoneproject).^p5).*((HH_3./sumHHzoneproject).^p6).*((HH_4./sumHHzoneproject).^p7).*((HH_5./sumHHzoneproject).^p8).*((HH_6./sumHHzoneproject).^p9).*COMMzoneproject;
                end

                
                Ldable= Lr_unrestricted+ Lb_unrestricted+ Lc_unrestricted+ landsP;
                %Enforce that Lr+Lb+Lc <= landddable
                Lrproject=Lr_unrestricted;
                
                Lbproject=Lb_unrestricted;
                
                Lcproject=Lc_unrestricted;

                L_change=zeros(dimZone,1);

                for i=1:dimZone;
                    if Ldable (i,1) > landdable (i,1);
                        Lrproject(i,1)=Lr_unrestricted(i,1).*(landdable(i,1)./Ldable(i,1));
                        Lbproject(i,1)=Lb_unrestricted(i,1).*(landdable(i,1)./Ldable(i,1));
                        Lcproject(i,1)=Lc_unrestricted(i,1).*(landdable(i,1)./Ldable(i,1));
                        L_change(i,1)=1;
                    end;
                end;

               % Adding the residuals ...
                if(count == 1)
                temp_Lr = Lrproject + eps_LU(:,1);
                temp_Lb = Lbproject + eps_LU(:,2);
                temp_Lc = Lcproject + eps_LU(:,3);
                end;
                if(count == 2)
                temp_Lr = Lrproject + (0.75*eps_LU(:,1));
                temp_Lb = Lbproject + (0.75*eps_LU(:,2));
                temp_Lc = Lcproject + (0.75*eps_LU(:,3));     
                    end;
                if(count == 3)
                temp_Lr = Lrproject + (0.5*eps_LU(:,1));
                temp_Lb = Lbproject + (0.5*eps_LU(:,2));
                temp_Lc = Lcproject + (0.5*eps_LU(:,3));  
                end;
                if(count == 4)
                temp_Lr = Lrproject + (0.25*eps_LU(:,1));
                temp_Lb = Lbproject + (0.25*eps_LU(:,2));
                temp_Lc = Lcproject + (0.25*eps_LU(:,3));        
                end;

                if(count == 5)
                temp_Lr = Lrproject;
                temp_Lb = Lbproject;
                temp_Lc = Lcproject;        
                end;

                for i=1:dimZone;
                    if(temp_Lr(i,1) > 0 ) Lrproject(i,1)= temp_Lr(i,1);end;
                    if(temp_Lb(i,1) > 0 ) Lbproject(i,1)= temp_Lb(i,1);end;
                    if(temp_Lc(i,1) > 0 ) Lcproject(i,1)= temp_Lc(i,1);end;
                    if(isnan(Lrproject(i,1)) == 1) Lrproject(i,1) = landr(i);end;
                    if(isnan(Lbproject(i,1)) == 1) Lbproject(i,1) = landb(i);end;
                    if(isnan(Lcproject(i,1)) == 1) Lcproject(i,1) = landc(i);end;
                end;        

%}
                for i=1:dimZone;
                    if(isnan(Lrproject(i,1)) == 1)   Lrproject(i,1) = Lrproject(i-1,1);end;
                    if(isnan(Lbproject(i,1)) == 1)   Lbproject(i,1) = Lbproject(i-1,1);end;
                    if(isnan(Lcproject(i,1)) == 1)   Lcproject(i,1) = Lcproject(i-1,1);end;
                        %temp_i = i
                       %temp_value = Lcproject(i,1)
                end;

%%%%%%%%%%%%%%%%%%%%% END OF LUDENSITY %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%% PRINTING THE VALUES %%%%%%%%%%%%%%%%%%%%%%
                if(count ==1 )

                csvwrite('C:\G-LUM Finale\Prediction\Output\EMP_timepd1.dat',EMPproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\HH_timepd1.dat',HHproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lr_timepd1.dat',Lrproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lb_timepd1.dat',Lbproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lc_timepd1.dat',Lcproject);

                end;


                if(count ==2)

                csvwrite('C:\G-LUM Finale\Prediction\Output\EMP_timepd2.dat',EMPproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\HH_timepd2.dat',HHproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lr_timepd2.dat',Lrproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lb_timepd2.dat',Lbproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lc_timepd2.dat',Lcproject);
                end;


                if(count ==3 )

                csvwrite('C:\G-LUM Finale\Prediction\Output\EMP_timepd3.dat',EMPproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\HH_timepd3.dat',HHproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lr_timepd3.dat',Lrproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lb_timepd3.dat',Lbproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lc_timepd3.dat',Lcproject);

                end;


                if(count ==4 )

                csvwrite('C:\G-LUM Finale\Prediction\Output\EMP_timepd4.dat',EMPproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\HH_timepd4.dat',HHproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lr_timepd4.dat',Lrproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lb_timepd4.dat',Lbproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lc_timepd4.dat',Lcproject);

                end;


                if(count ==5 )
                csvwrite('C:\G-LUM Finale\Prediction\Output\EMP_timepd5.dat',EMPproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\HH_timepd5.dat',HHproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lr_timepd5.dat',Lrproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lb_timepd5.dat',Lbproject);
                csvwrite('C:\G-LUM Finale\Prediction\Output\Lc_timepd.dat',Lcproject);

                end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Jen's TDM and PRINTING VALUES (Tried to Integrate with Brenda's code by Siva)%%%%%%%%%%%%%%%%%

%        EMP_all(dimZone*(dimTimePd*(jj-1)+count-1)+1:dimZone*(dimTimePd*(jj-1)+count),:)=EMPproject;        
%        HH_all(dimZone*(dimTimePd*(jj-1)+count-1)+1:dimZone*(dimTimePd*(jj-1)+count),:)=HHproject;  
%        Lr_all(dimZone*(dimTimePd*(jj-1)+count-1)+1:dimZone*(dimTimePd*(jj-1)+count),:)=Lrproject;  
%        Lb_all(dimZone*(dimTimePd*(jj-1)+count-1)+1:dimZone*(dimTimePd*(jj-1)+count),:)=Lbproject;  
%        Lc_all(dimZone*(dimTimePd*(jj-1)+count-1)+1:dimZone*(dimTimePd*(jj-1)+count),:)=Lcproject;  
%        c1 = fourstep(Lrproject,Lbproject,Lcproject,HHproject,EMPproject,c,X,Y,Z,count);
%        c=c1;
%        %c_all(dimLink*(dimTimePd*(jj-1)+count-1)+1:dimLink*(dimTimePd*(jj-1)+count),:) = c;
%        vc = load('C:\G-LUM Finale\Prediction\Output\V-over-C.dat'); 
%        vc_all(dimLink*(dimTimePd*(jj-1)+count-1)+1:dimLink*(dimTimePd*(jj-1)+count),:) = vc;
%        
%        if (count == dimTimePd)
%            link_speed = load('C:\G-LUM Finale\Prediction\Output\link_speed.dat'); 
%            link_length = load('C:\G-LUM Finale\Prediction\Input\link_length.txt');
%            weight = link_length;
%            mean_speed = (sum(weight.*link_speed))/sum(weight);
%            variance_speed = (sum(weight.*link_speed.*link_speed)*sum(weight) - (sum(weight.*link_speed))^2)/((sum(weight))^2-sum(weight.*weight));
%            speed_all(jj,:) = horzcat(mean_speed,variance_speed);
%        end;
        

%%%%%%%%%%%%%%%% REINIATIALIZATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Assign the predicted values as base year values to continue ...
                    EMPbase_original = EMPproject;
                    HHbase_original = HHproject;
                    landr = Lrproject;
                    landb = Lbproject;
                    landc = Lcproject;

                    landd = landr + landb + landc + landsP; % The streets is assumed to remain the same ..
                    landdable = ZoneSize - landeP; % Developable land
                    landv = ZoneSize - landd - landeP;
                for i=1:dimZone
                    if(landr(i) < 0.0001) landr(i) = 0.0001;end;
                    if(landb(i) < 0.0001) landb(i) = 0.0001;end;
                    if(landc(i) < 0.0001) landc(i) = 0.0001;end;
                    if(landv(i) < 0.0001) landv(i) = 0.0001;end;
                    if (landdable(i)<0.001) landdable(i)=0.001;end; 
                end;

landdRatio = (landr + landb  + landc + landsP )./( landdable);  % Proportion of developable land developed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                for i=1:dimZone;
                    for j = 1:dimEMP
                        if(isreal(EMPproject(i,j) ) == 0 )
                            temp_EMP = j
                        end;
                    end;

                    for k = 1:dimHH
                        if(isreal(HHproject(i,k))==0) 
                            temp_HHi = i 
                        end;
                    end;
                end;

        end; % end of count 
    landrPproject = Lrproject;
    landbPproject = Lbproject;
    landcPproject= Lcproject;
    landsPproject = landsP;
    landePproject = landeP;
    landuPproject = ZoneSize - Lrproject - Lcproject - Lbproject - landsP - landeP; 

    landall=horzcat(landrPproject,landbPproject,landePproject,landuPproject,landcPproject,landsPproject);
end; % ii loop 

%csvwrite('C:\G-LUM Finale\Prediction\Output\EMP_all.dat',EMP_all);
%csvwrite('C:\G-LUM Finale\Prediction\Output\HH_all.dat',HH_all);
%csvwrite('C:\G-LUM Finale\Prediction\Output\Lr_all.dat',Lr_all);
%csvwrite('C:\G-LUM Finale\Prediction\Output\Lb_all.dat',Lb_all);
%csvwrite('C:\G-LUM Finale\Prediction\Output\Lc_all.dat',Lc_all);
%csvwrite('C:\G-LUM Finale\Prediction\Output\vc_all.dat',vc_all);
%csvwrite('C:\G-LUM Finale\Prediction\Output\speed_all.dat',speed_all);
sprintf('Finished Predicting from base year to %d time period(s) in future.',dimTimePd)


function edit1_Callback(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of edit1 as text
%        str2double(get(hObject,'String')) returns contents of edit1 as a double
years = str2double(get(hObject,'string'));
if isnan(years)
	errordlg('You must enter a numeric value','Bad Input','modal')
	return
end
csvwrite('C:\G-LUM Finale\Prediction\Output\years.txt',years);

% --- Executes during object creation, after setting all properties.
function edit1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end  

