%move simple test
%note: interesting things happen if alpha is large (.75)
work=[];
penss=[];
ww=[];
hphp=[];
relax=.1;
setvalues4
sumL=sum(L);
sumpens=sum(pens);
theta_work=zeros(N,1);
theta_pens=zeros(N,1);
A=ones(N,1);
agl=0;
%(1) symmetry
tau=[1 1.1 1.1;1.1 1 1.1; 1.1 1.1 1]
oldw=ones(N,1);
oldact_work=ones(N,1);
oldact_pens=ones(N,1);
goon=1;
while (goon)
 expend;
 [optw,msg]=fsolve("xssw",w(2:N));
 w=[1;optw];
 expend;
 Hp=house(H);
 p=price();
 [util_work,util_pens]=util(p,Hp);
 loc_work=locshare(util_work,theta_work);
 loc_pens=locshare(util_pens,theta_pens);
 act_work=L/sum(L);
 act_pens=pens/sum(pens);
 norm=(w-oldw)'*(w-oldw);
 norm=max(norm,(oldact_work-act_work)'*(oldact_work-act_work));
 norm=max(norm,(oldact_pens-act_pens)'*(oldact_pens-act_pens));
 if (norm<1e-20)
   goon=0;
 end
 L=sumL*(relax*loc_work+(1-relax)*act_work);
 pens=sumpens*(relax*loc_pens+(1-relax)*act_pens);
 pensshare=pens/sum(pens);
 oldw=w;
 oldact_work=act_work;
 oldact_pens=act_pens;
end
act_work
act_pens
kshare=K/sum(K)
w
rentincome=sum(Hp.*H)+r*sumK;
workerincome=w+(1-eta)*rentincome/sum(L)
pensincome=eta*rentincome/sum(pens)
check=L.*workerincome+pens.*pensincome-expenditures
shar=share()
%(2) isolated region
tau=[1 1.1 1.3;1.1 1 1.3; 1.3 1.3 1]
oldw=ones(N,1);
oldact_work=ones(N,1);
oldact_pens=ones(N,1);
goon=1;
while (goon)
 expend;
 [optw,msg]=fsolve("xssw",w(2:N));
 w=[1;optw];
 expend;
 Hp=house(H);
 p=price();
 [util_work,util_pens]=util(p,Hp);
 loc_work=locshare(util_work,theta_work);
 loc_pens=locshare(util_pens,theta_pens);
 act_work=L/sum(L);
 act_pens=pens/sum(pens);
 norm=(w-oldw)'*(w-oldw);
 norm=max(norm,(oldact_work-act_work)'*(oldact_work-act_work));
 norm=max(norm,(oldact_pens-act_pens)'*(oldact_pens-act_pens));
 if (norm<1e-20)
   goon=0;
 end
 L=sumL*(relax*loc_work+(1-relax)*act_work);
 pens=sumpens*(relax*loc_pens+(1-relax)*act_pens);
 pensshare=pens/sum(pens);
 oldw=w;
 oldact_work=act_work;
 oldact_pens=act_pens;
end
act_work
act_pens
kshare=K/sum(K)
w
rentincome=sum(Hp.*H)+r*sumK;
workerincome=w+(1-eta)*rentincome/sum(L)
pensincome=eta*rentincome/sum(pens)
check=L.*workerincome+pens.*pensincome-expenditures
shar=share()
%(3) increase pensioners
inhab=L+pens;
printf("inhab L pens w p Hp uwork upens");
[inhab L pens w p Hp util_work util_pens] 
sumpens=sumpens*1.40;
sumL=sumL*.8;
eta=eta*1.4;
goon=1;
oldw=ones(N,1);
oldact_work=ones(N,1);
oldact_pens=ones(N,1);
while (goon)
 expend;
 [optw,msg]=fsolve("xssw",w(2:N));
 w=[1;optw];
 expend;
 Hp=house(H);
 p=price();
 [util_work,util_pens]=util(p,Hp);
 loc_work=locshare(util_work,theta_work);
 loc_pens=locshare(util_pens,theta_pens);
 act_work=L/sum(L);
 act_pens=pens/sum(pens);
 norm=(w-oldw)'*(w-oldw);
 norm=max(norm,(oldact_work-act_work)'*(oldact_work-act_work));
 norm=max(norm,(oldact_pens-act_pens)'*(oldact_pens-act_pens));
 if (norm<1e-20)
   goon=0;
 end
 L=sumL*(relax*loc_work+(1-relax)*act_work);
 pens=sumpens*(relax*loc_pens+(1-relax)*act_pens);
 pensshare=pens/sum(pens);
 oldw=w;
 oldact_work=act_work;
 oldact_pens=act_pens;
end
inhab=L+pens;
[inhab L pens w p Hp util_work util_pens]
act_work
act_pens
kshare=K/sum(K)
w
rentincome=sum(Hp.*H)+r*sumK;
workerincome=w+(1-eta)*rentincome/sum(L)
pensincome=eta*rentincome/sum(pens)
check=L.*workerincome+pens.*pensincome-expenditures
shar=share()