%2018.10.01: This is to illustrate the process to apportion the number of families for 9 building classes into 17 building sub-types.
clear;
clc;
load isomat_V2.mat
isomat = isomat_V2;
strucout=isomat(:,8:16);
strucout2=isomat(:,8:16); % number of families sorted by: 1-storey,2-3,4-6,7-9,10,steel/RC,mixed masonry,brick/wood,other
line_no = length(isomat(:,1));
pop_family = isomat(:,4); % average population per family
building_sub_types = 17;
% % 'BRIWOMC1','BRIWOMC23','STLRCMC1','STLRCMC23','STLRCMC46','STLRCMC79','STLRCMC10', ...
% % 'MIXEDMC1','MIXEDMC23','MIXEDMC46','MIXEDMC79','MIXEDMC10','OTHERMC1','OTHERMC23','OTHERMC46','OTHERMC79','OTHERMC10'
% resmatrix = zeros(line_no,17);
% C1-C9 corresponds to famailies living in 1-storey,2-3,4-6,7-9,G10,steel/RC,mixed masonry,brick/wood,other
%% brick/wood (C8) and 1-story (C1)
tmp1 = min(strucout2(:,1),strucout(:,8)); % should be equal to resmatrix(:,1)

resmatrix(:,1)=strucout(:,8)-strucout(:,1); 
replacewood1=find(resmatrix(:,1)<0); % C8 < C1 (negative)
replacewood1a=find(resmatrix(:,1)>=0); % C8 >= C1 (positive)

strucout(replacewood1,1)=strucout(replacewood1,1)-strucout(replacewood1,8); % C1-C8
resmatrix(replacewood1,1)=strucout(replacewood1,8); % min(C1,C8) = C8
strucout(replacewood1,8)=0;
%when brick/wood buildings are more than 1-storey buildings
strucout(replacewood1a,8)=strucout(replacewood1a,8)-strucout(replacewood1a,1);% C8-C1
strucout(replacewood1a,1)=0;
resmatrix(replacewood1a,1)=strucout2(replacewood1a,1); % min(C1,C8) = C1
%Final result
%strucout(:,1): C1-C8 & 0
%strucout(:,8): 0 & C8-C1 
%resmatrix(:,1): min(C1,C8)

%% brick/wood (C8 -> "C8-C1 or 0") and 2-3 storey (C2)
tmp2 = min(strucout2(:,2),strucout(:,8));% should be equal to resmatrix(:,2)

resmatrix(:,2)=strucout(:,8)-strucout(:,2);
replacewood2=find(resmatrix(:,2)<0); % C8-C1 < C2 & 0 < C2
replacewood2a=find(resmatrix(:,2)>=0); % C8-C1 > C2
strucout(replacewood2,2)=strucout(replacewood2,2)-strucout(replacewood2,8); % C2 - (C8-C1) or C2 - 0
resmatrix(replacewood2,2)=strucout(replacewood2,8); % min(C2,C8-C1)=C8-C1 or min(C2,0)=0
strucout(replacewood2,8)=0; % 
%when remaining brick/wood is greater than 2-3 storey
strucout(replacewood2a,8)=strucout(replacewood2a,8)-strucout(replacewood2a,2);% (C8-C1) - C2
strucout(replacewood2a,2)=0;
resmatrix(replacewood2a,2)=strucout2(replacewood2a,2); % min(C8-C1,C2) = C2
%Final result
%strucout(:,2): C2-(C8-C1) or C2 & 0
%strucout(:,8): 0 & (C8-C1)-C2
% resmatrix(:,2): min(C2,C8-C1) or min(C2,0)=0

%% steel/RC (C6) and G10 storey (C5)
tmp7 = min(strucout2(:,5),strucout(:,6)); % should be equal to resmatrix(:,7)

resmatrix(:,7)=strucout(:,6)-strucout(:,5);
replacerc10=find(resmatrix(:,7)<0); % C6 < C5
replacerc10a=find(resmatrix(:,7)>=0); % C6 >= C5
strucout(replacerc10,5)=strucout(replacerc10,5)-strucout(replacerc10,6); % C5-C6
resmatrix(replacerc10,7)=strucout(replacerc10,6); % min(C5,C6) = C6
strucout(replacerc10,6)=0;
%when steel/RC is greater than G10 storey
strucout(replacerc10a,6)=strucout(replacerc10a,6)-strucout(replacerc10a,5); % C6-C5
strucout(replacerc10a,5)=0;
resmatrix(replacerc10a,7)=strucout2(replacerc10a,5); % min(C5,C6) = C5
%Final result
%strucout(:,5): C5-C6 & 0
%strucout(:,6): 0 & C6-C5
% resmatrix(:,7): min(C5,C6)

%% steel/RC (C6 -> "0 or C6-C5") and 7-9 storey (C4) 
tmp6 = min(strucout2(:,4),strucout(:,6)); % should be equal to resmatrix(:,6)

resmatrix(:,6)=strucout(:,6)-strucout(:,4);
replacerc9=find(resmatrix(:,6)<0); % 0 < C4 or C6-C5 < C4
replacerc9a=find(resmatrix(:,6)>=0); % C6-C5 >= C4
strucout(replacerc9,4)=strucout(replacerc9,4)-strucout(replacerc9,6); % C4-0 or C4-(C6-C5)
resmatrix(replacerc9,6)=strucout(replacerc9,6); % min(0,C4)=0 or min(C6-C5,C4) = C6-C5
strucout(replacerc9,6)=0;
%when remaining steel/RC is greater than 7-9 storey
strucout(replacerc9a,6)=strucout(replacerc9a,6)-strucout(replacerc9a,4); % (C6-C5)-C4
strucout(replacerc9a,4)=0;
resmatrix(replacerc9a,6)=strucout2(replacerc9a,4); % min(C6-C5,C4) = C4
%Final result
%strucout(:,4): C4 or C4-(C6-C5) & 0
%strucout(:,6): 0 & (C6-C5)-C4
% resmatrix(:,6): min(0,C4)=0 or min(C6-C5,C4)

%% steel/RC (C6 -> "0 or (C6-C5)-C4") and 4-6 storey (C3) 
tmp5 = min(strucout2(:,3),strucout(:,6)); % should be equal to resmatrix(:,5)

resmatrix(:,5)=strucout(:,6)-strucout(:,3);
replacerc46=find(resmatrix(:,5)<0); % 0<C3 or (C6-C5)-C4 < C3
replacerc46a=find(resmatrix(:,5)>=0); % (C6-C5)-C4 >= C3
strucout(replacerc46,3)=strucout(replacerc46,3)-strucout(replacerc46,6); % C3-0 or C3-[(C6-C5)-C4]
resmatrix(replacerc46,5)=strucout(replacerc46,6); % min(0,C3)=0 or min([(C6-C5)-C4],C3) = (C6-C5)-C4
strucout(replacerc46,6)=0;
%when remaining steel/RC is greater than 4-6 storey
strucout(replacerc46a,6)=strucout(replacerc46a,6)-strucout(replacerc46a,3); % [(C6-C5)-C4]-C3
strucout(replacerc46a,3)=0;
resmatrix(replacerc46a,5)=strucout2(replacerc46a,3); % min([(C6-C5)-C4],C3) = C3
%Final result
%strucout(:,3): C3 or C3-[(C6-C5)-C4] & 0
%strucout(:,6): 0 & [(C6-C5)-C4]-C3
% resmatrix(:,5): min(0,C3)=0 or min([(C6-C5)-C4],C3)

%% steel/RC (C6 -> "0 or [(C6-C5)-C4]-C3")  and 2-3 storey (C2 -> "C2-(C8-C1) or C2 or 0")
tmp4 = min(strucout(:,2),strucout(:,6)); % should be equal to resmatrix(:,4)

resmatrix(:,4)=strucout(:,6)-strucout(:,2);
replacerc23=find(resmatrix(:,4)<0);
replacerc23a=find(resmatrix(:,4)>=0); %
strucout(replacerc23,2)=strucout(replacerc23,2)-strucout(replacerc23,6);
resmatrix(replacerc23,4)=strucout(replacerc23,6); 
strucout(replacerc23,6)=0;
%when remaining steel/RC is greater than 2-3 storey
strucout(replacerc23a,6)=strucout(replacerc23a,6)-strucout(replacerc23a,2); % {[(C6-C5)-C4]-C3} - [C2-(C8-C1)] or {[(C6-C5)-C4]-C3} - C2 or {[(C6-C5)-C4]-C3}
strucout(replacerc23a,2)=0;
resmatrix(replacerc23a,4)=strucout(replacerc23a,2); %
%% steel/RC (C6 -> "{[(C6-C5)-C4]-C3} - [C2-(C8-C1)] or {[(C6-C5)-C4]-C3} - C2 or {[(C6-C5)-C4]-C3} or 0") and 1-storey (C1 -> "C1-C8 or 0") 
tmp3 = min(strucout(:,1),strucout(:,6));% should be equal to resmatrix(:,3)

resmatrix(:,3)=strucout(:,6)-strucout(:,1);
replacerc1=find(resmatrix(:,3)<0);
replacerc1a=find(resmatrix(:,3)>=0); %Case 2 brick/wood is larger than 1-storey
strucout(replacerc1,1)=strucout(replacerc1,1)-strucout(replacerc1,6);
resmatrix(replacerc1,3)=strucout(replacerc1,6); 
strucout(replacerc1,6)=0;
%when remaining steel/RC is greater than 2-3 storey
strucout(replacerc1a,6)=strucout(replacerc1a,6)-strucout(replacerc1a,1);
strucout(replacerc1a,1)=0;
resmatrix(replacerc1a,3)=strucout(replacerc1a,1); %

% Assign "mixed masonry" and "other" buildings into 1/2-3/4-6/7-9/G10 storey classes 
for j=1:5
    % "mixed masonry":
    resmatrix(:,7+j) = strucout(:,j).*( strucout(:,7)./(strucout(:,7)+strucout(:,9)) );
    % "other":
    resmatrix(:,12+j) = strucout(:,j).*( strucout(:,9)./(strucout(:,7)+strucout(:,9)) );
end
% resmatrix(:,8:12)=strucout(:,1:5).* (strucout(:,7)./(strucout(:,7)+strucout(:,9)) );
% resmatrix(:,13:17)=strucout(:,1:5).* (strucout(:,9)./(strucout(:,7)+strucout(:,9)) );

resmatrix(isnan(resmatrix))=0;
check2 = zeros(93,9);

%% The following part is to compare the apportioned buildings with census records.

% The columns of "check2" represent the derived number of families (from discretization results) living in
% each building class; these can be compared with the corresponding columns
% in "strucout2"
check2(:,1) = resmatrix(:,1)+resmatrix(:,3)+resmatrix(:,8)+resmatrix(:,13); % 1-storey
check2(:,2) = resmatrix(:,2)+resmatrix(:,4)+resmatrix(:,9)+resmatrix(:,14); % 2-3 storey
check2(:,3) = resmatrix(:,5)+resmatrix(:,10)+resmatrix(:,15); % 4-6 storey
check2(:,4) = resmatrix(:,6)+resmatrix(:,11)+resmatrix(:,16); % 7-9 storey
check2(:,5) = resmatrix(:,7)+resmatrix(:,12)+resmatrix(:,17); % â‰?10 storey (G10)
check2(:,6) = resmatrix(:,3)+resmatrix(:,4)+resmatrix(:,5)+resmatrix(:,6)+resmatrix(:,7); % "steel/RC"
check2(:,7) = resmatrix(:,8)+resmatrix(:,9)+resmatrix(:,10)+resmatrix(:,11)+resmatrix(:,12); % "mixed masonry"
check2(:,8) = resmatrix(:,1)+resmatrix(:,2); % "brick/wood"
check2(:,9) = resmatrix(:,13)+resmatrix(:,14)+resmatrix(:,15)+resmatrix(:,16)+resmatrix(:,17); % "other"

sum(check2,1) % sum of apportioned number of families
sum(strucout2,1) % census recorded number of families

% % original_pop_in_house = zeros(93,9); % originally recorded population in 2010 census
% % check_pop_in_house = zeros(93,9); % population derived from the 17 sub-types9 buildings, as apportioned from the 9 building classes above
% % for k=1:93
% %     check_pop_in_house(k,:) = check2(k,:)*10*isomat(k,17)*pop_family(k,1);
% %     original_pop_in_house(k,:) = strucout2(k,:)*10*isomat(k,17)*pop_family(k,1); 
% % end
% % sum(check_pop_in_house,1) % sum of apportioned population
% % sum(original_pop_in_house,1) % census recorded population