% By Paul Dreik 20210108
% License GPL v3
% https://www.pauldreik.se/

clear,clc, close all hidden

%product, manufacturer, time in minutes, selling price on trade depot
products={
'doughnut','donut shop',45,950';

'grass','gardening supplies',30,310;
'tree saplings','gardening supplies',90,420;

'hammer','hardware store',14,90;
'measuring tape','hardware store',20,110;
'shovel','hardware store',30,150;
'cooking utensils','hardware store',45,250;

'chair','furniture store',20,300;
'table','furniture store',30,500;

'nails','building supplies store',5,80;
'planks','building supplies store',30,120;
'bricks','building supplies store',20,0;
'cement','building supplies store',50,440;
'glue','building supplies store',60,440;

'vegetables','farmers market',20,160;
'flour bag','farmers market',30,570;
'fruit and berries','farmers market',90,730;

'wood','factory',3,0;
'metal','factory',1,0;
'sugar and spices','factory',240,110;
'textiles','factory',180,90;
'plastic','factory',9,0;
'seeds','factory',20,0;
'minerals','factory',30,0;
'chemicals','factory',120,60};


%how many items can be produced per hour, assuming infinite
% supply of dependencies
fraction_per_hour=60./cell2mat(products(:,3));
prices=cell2mat(products(:,4));

factories=unique(products(:,2));

% element i,j says how many items of type j are needed
% to produce an item i. I is the adjacency matrix
% of the graph.
Deps=[
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	1	0	0	0	0	0
0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0
0	0	0	0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	0	1	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	1	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	1	0	0	1	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	2	0	0	2	0	0	0
0	0	0	1	0	0	0	0	0	1	0	0	0	0	0	0	0	2	0	0	0	0	0	0	0
0	0	0	1	0	0	0	0	0	2	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	1
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	0	0	2
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	2	0	0
0	0	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	2	0	0
0	0 0 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0 0 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0 0 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0 0 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0 0 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0 0 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0 0 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
0	0 0 0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
];
% the self consumed production is (production'*Deps)

effort_choice=[
1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0 % doughnut shop
0	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0 % gardening
0	0	0	1	1	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0 % hardware
0	0	0	0	0	0	0	1	1	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0 % furniture
0	0	0	0	0	0	0	0	0	1	1	1	1	1	0	0	0	0	0	0	0	0	0	0	0 % build supply
0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	1	1	0	0	0	0	0	0	0	0 % farmer
0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	1	1	1	1	1	1	1	1 % factory
];

Nproducts=size(Deps,1);

% required factory capacity is effort_matrix*production
effort_matrix=effort_choice*diag(1./fraction_per_hour);
Nfactories=size(effort_matrix,1);
% how high may the effort be?
% we have 22 parallel tracks total in the basic factory, all others (doughnut etc) are single track.
max_effort=ones(Nfactories,1);
max_effort(end)=22;

%production is the number of items per hour produced of each kind.
%net_production is what remains after some of the production is used for producing
% other products
%
% maximize profit=prices'*net_production
% where net_production=production-production'*Deps
% subject to not exceeding capacity of any factory:
% B*production<=1for all elements
%
% let x=production, all elements >=0.
% maximize: prices'*x*(I-Deps) = ((I-Deps)*prices)'*x
% with constraints B*x <=1
%
ID=eye(size(Deps))-Deps;
c = ID*prices;
A = [effort_matrix;
ID'];
b = [max_effort;
zeros(Nproducts,1)];
lb = zeros(Nproducts,1);
ub = ones(Nproducts,1)*10000000000000;
ctype = [repmat('U',1,Nfactories) repmat('L',1,Nproducts)]; % upper bound. lower bound
vartype =repmat('C',Nproducts,1); % continuous
s = -1; % maximization

param.msglev = 1;
param.itlim = 100;

[xmin, fmin, status, extra] = ...
glpk (c, A, b, lb, ub, ctype, vartype, s, param);
netprod=ID'*xmin;
for i=1:Nproducts
  printf('produce %.3f of %s, net prod %.3f\n',xmin(i),products{i,1},netprod(i));     
endfor

fid=fopen('results.table.md','w');
assert(fid>0);
fprintf(fid,'|  Item  |   Production [items/hour] |   Net production [items/hour] | Profit [coins/hour] |\n');
fprintf(fid,'|-----|-----|-----|-----|\n');
for i=1:Nproducts
a=  xmin(i);
b=netprod(i);
if abs(a)<1e-3
  a=0;  
endif
if abs(b)<1e-3
  b=0;
endif
fprintf(fid,'| %s | %g | %g | %g |\n',products{i,1},a,b,b*prices(i));     
endfor
fclose(fid);

%export to a graphviz graph
fid=fopen('simcity.dot','w');
assert(fid>0);
fprintf(fid,'digraph G {\n');
for i=1:length(factories)  
  fprintf(fid,'subgraph cluster_%s {\n',strrep(factories{i},' ','_'));
  fprintf(fid,'label="%s"\n',factories{i});
  fprintf(fid,'color=blue\n');
  I=find(ismember(products(:,2),factories{i}));
  for j=I(:)'
    fprintf(fid,'%s [label="%s"]\n',strrep(products{j,1},' ','_'),products{j,1});
  endfor
  fprintf(fid,'}\n');
endfor
for i=1:size(Deps,1)
  for j=1:size(Deps,1)
    if Deps(i,j)
      fprintf(fid,'%s -> %s\n',strrep(products{j,1},' ','_'),strrep(products{i,1},' ','_'));
    endif
  endfor 
endfor
fprintf(fid,'}\n');

fclose(fid);