function [gurobi_model,linmodel] = srec_formulate_mip_new( samplegraph,mipopt )
% formulate the whole problem into a mixed integer programing
% variables:
% x=[x1,x2...xn,    %node selection, xk=1 iff kth node is selected
%   y1,y2,...ym,    %edge selection, yk=1 iff edge tril(adjmat)(k) is selected
%   f1,f2,...f_2m,  %flow assignment, used to ensure connectivity
%   g1,g2,...gn,    %flow source, used to ensure connectivity
%   u1,u2,...us,    %2d segment coverage, uk=1 iff kth segment is
%                   hard-assigned to some element
%   c1,c2,...cs,    %2d segment overlap count, ck>0 iff kth segment is 
%                   hard-assigned ck+1 times
%   r11,r12...r1t,r21...r_nk,  %2d segment overlap for each component,
%                               r_ij=1 if jth 2d segment for ith component
%                               is overlapped by others (hard-assigned to
%                               others)
%   h11,h12...hmn,  %label assignment to 2d segment, hij=1 iff ith segment 
%                    is assigned jth label in dblbs (hard-assigned to a
%                    component of label[j])
%   t11,t12,...tij...,    %group label activation indicator, tij = 1 if
%                          at least one of the 2d segments in ith group take jth label
%   v11,v12,...vnk, %2d segment hard assignment, vij=1 iff ith segment IS
%                    assigned to jth component (i.e. covered by jth
%                    component)

% weights
w_cover = mipopt.m_weight_cover; %best = 100
w_mustcover = mipopt.m_weight_mustcover; %best = 50
w_match = mipopt.m_weight_match; %best = 10
w_overlap = mipopt.m_weight_overlap; %best = 15
w_struct = mipopt.m_weight_struct; %best = 2
w_segroup = mipopt.m_weight_segroup; %best = 5
w_hardassign = mipopt.m_weight_hardassign;

% quantize match and struct score?
% set to 0 to disable quantization
nbin_match = mipopt.m_nbin_unary;
nbin_struct = mipopt.m_nbin_pairwise;

% global forbiden elements
%{
[~,objtype] = util_get_name_by_label(samplegraph.m_dblbs(1));
objtype = objtype{1};
segids_forbid_always = srec_forbid_special_component(objtype);
%}
segids_forbid_always = [];

% if isequal(objtype,'chair')
%     segids_forbid_always = [13 3;18 1]';
% else
%     segids_forbid_always = [37 3]';
% end

% for each selected element, the coverlen (proportion covered) must be larger than this
thres_cover_ratio_sufficient = 0.8;

% for each selected element, if the coverlen is below this threshold,
% generates a penalty, otherwise generates an award
thres_min_coverlen_to_award = 1;

max_overlap_ratio = 0.9;

if nargin < 2
    mipopt = SBSRMIPOPT;
end

use_overlap_hard_cons = mipopt.m_useoverlap;

% for each segment, which component covers it?
nseg = max(samplegraph.m_seglbs_per_sample);
ncomp = length(samplegraph.m_complabels);
seglbs2comp = false(nseg,ncomp);

switch mipopt.m_covermode
    case 'normal'
        comp2sample = samplegraph.m_adjmat_comp2sample;
    case 'fullcover'
        comp2sample = samplegraph.m_adjmat_comp2sample_fullcover;
    case 'justenough'
        comp2sample = samplegraph.m_adjmat_comp2sample_justenough;
end

for i=1:ncomp
    seglbs = samplegraph.m_seglbs_per_sample(comp2sample(i,:));
    seglbs2comp(seglbs,i) = true;
end
idxmap_seglbs2comp = double(seglbs2comp); %give each seg2comp a linear index
idxmap_seglbs2comp(seglbs2comp) = 1:nnz(seglbs2comp);

% ignore f,g first, the variable is x=[x1..xn,y1...ym,u1...us]
if mipopt.m_weight_struct == 0
    adjmat = samplegraph.m_structscore_comp2comp_raw > 0;
else
    % note, scored graph no longer filters invalid edges
    % you must filter it explicitly
    adjmat = samplegraph.get_scored_graph>0;
    adjmat = adjmat & ~samplegraph.m_adjmat_invalid_comp2comp;
    
%     if mipopt.m_useoverlap
%         adjmat = samplegraph.get_scored_graph>0; %not filtering overlap edges
%     else
%         % overlap edges explicitly filtered
%         adjmat = samplegraph.get_scored_graph('adjmat',~samplegraph.m_adjmat_invalid_comp2comp)>0;
%     end
end
adjmat = tril(adjmat);
nx = size(adjmat,1);
ny = nnz(adjmat);
nu = max(samplegraph.m_seglbs_per_sample);
nc = nu;
nr = nnz(seglbs2comp);
if ~use_overlap_hard_cons
    nr = 0;
end
nv = nnz(seglbs2comp);

%% segment continuous grouping
segcontgroup = mipopt.m_segcontgroup;
contgroup2seglbs = false(length(segcontgroup),size(seglbs2comp,1)); %m(i,j)==1 iff segment j belongs to continuous group i
for i=1:length(segcontgroup)
    seglbs = segcontgroup{i};
    contgroup2seglbs(i,seglbs) = true;
end

%% used to introduce segment group cost
nh = 0; %hij=1 iff ith segment takes jth label
nt = 0; %tij=1 iff some segment in ith group takes jth label
seglbsgroup = {}; %seglbsgroup{i} = seglbs for ith group

% label_per_seglbsgroup{i} = [label1,label2..labelk], if one segment in ith
% group takes a label in this set, other segments preferrably take the
% same label
label_per_seglbsgroup = {};
hidx_seglbs2label = [];
    
if ~isempty(mipopt.m_seglbsgroup)
    if any(size(mipopt.m_seglbsgroup)==1) %no group label
        ngroup = length(mipopt.m_seglbsgroup);
        seglbsgroup = mipopt.m_seglbsgroup;
        label_per_seglbsgroup = cell(1,ngroup);
        label_per_seglbsgroup(:) = {samplegraph.m_dblbs};
    else %have group label
        seglbsgroup = mipopt.m_seglbsgroup(1,:);
        label_per_seglbsgroup = mipopt.m_seglbsgroup(2,:);
    end
    
    % convert h's linear index into (seglb,label), that is,
    % h(i)->(seglb,label), meansh h(i) is for marking seglb with label
    maxseglb = size(seglbs2comp,1);
    [xx,yy] = ndgrid(1:maxseglb,samplegraph.m_dblbs);
    h2seglabelpair = cat(1,xx(:)',yy(:)');
    nh = size(h2seglabelpair,2);

    % map (seglb,label) to linear index of h
    hidx_seglbs2label =zeros(maxseglb,max(samplegraph.m_dblbs));
    idx = sub2ind(size(hidx_seglbs2label),h2seglabelpair(1,:),h2seglabelpair(2,:));
    hidx_seglbs2label(idx) = 1:length(idx);
end

% convert t's linear index into (groupid,label), that is,
% t(i)->(groupid,label) means t(i) is marking whether any of the
% segment in the group takes the label
t2grouplabelpair = zeros(2,0);
for i=1:length(label_per_seglbsgroup)
    thislbs = label_per_seglbsgroup{i};
    gpair = zeros(2,length(thislbs));
    gpair(1,:) = i;
    gpair(2,:) = thislbs;
    t2grouplabelpair = cat(2,t2grouplabelpair,gpair);
end
nt = size(t2grouplabelpair,2);

% map (groupid,label) to linear index of t
tidx_groupid2label = zeros(length(label_per_seglbsgroup),max(samplegraph.m_dblbs));
idx = sub2ind(size(tidx_groupid2label),t2grouplabelpair(1,:),t2grouplabelpair(2,:));
tidx_groupid2label(idx) = 1:length(idx);

%%
nvar = nx+ny+nu+nc+nr+nh+nt+nv;

% convert from variable to its index
x2varidx = 1:nx;
y2varidx = nx+(1:ny);
u2varidx = nx+ny+(1:nu);
c2varidx = nx+ny+nu+(1:nc);
r2varidx = nx+ny+nu+nc+(1:nr);
h2varidx = nx+ny+nu+nc+nr+(1:nh);
t2varidx = nx+ny+nu+nc+nr+nh+(1:nt);
v2varidx = nx+ny+nu+nc+nr+nh+nt+(1:nv);
x2rs = @(x)x2rs_(x,idxmap_seglbs2comp);
x2us = @(x)x2us_(x,seglbs2comp);
seglbs2vs = @(s)idxmap_seglbs2comp(s,idxmap_seglbs2comp(s,:)>0);

% if fixvars(i) is not nan, variable i is considered fixed
fixvars = nan(1,nvar);

% get vs that refers to certain sketch segment that is covered by elements
% of a given label
seglbs_label2vs = @(s,l)idxmap_seglbs2comp(s,...
    samplegraph.m_complabels==l & idxmap_seglbs2comp(s,:)>0);

% structure cost of each edge
structscore_per_edge = samplegraph.m_structscore_comp2comp_raw(adjmat);

if GLOBALCONST.DBG_CHECK('connect if cover','MIP: using inverse distance as structure score')
    distmat = squareform( pdist(samplegraph.m_compts'));
    structscore_per_edge = exp(mat2gray(distmat(adjmat)));
end

if ~isempty(structscore_per_edge)
    structscore_per_edge = mat2gray(structscore_per_edge);
    if nbin_struct
        structscore_per_edge = util_quantize(structscore_per_edge,nbin_struct);
    end
end
structcost_per_edge = 1-structscore_per_edge;
y2structcost = structcost_per_edge;
assert(length(y2structcost) == ny);

% 2d segment length
seglbs = unique(samplegraph.m_seglbs_per_sample);
lens = histc(samplegraph.m_seglbs_per_sample,seglbs);
lens_per_2dseg = zeros(1,nu);
lens_per_2dseg(seglbs) = lens;
weight_per_2dseg = lens_per_2dseg/mean(lens_per_2dseg);
weight_per_2dseg = weight_per_2dseg/max(weight_per_2dseg);
u2seglens = lens_per_2dseg;
u2segweight = weight_per_2dseg;
nsample = length(samplegraph.m_seglbs_per_sample);
u2ovcost = u2seglens/nsample;
c2ovcost = u2ovcost;

if ~isempty(mipopt.m_seglbs_use)
    % we have specified usable seglbs, other seglbs are considered covered
    % by ghost component
    seglbsuse = mipopt.m_seglbs_use;
    seglbsall = 1:max(samplegraph.m_seglbs_per_sample);
    seglbs_always_covered = setdiff(seglbsall,seglbsuse);
else
    seglbs_always_covered = [];
end

%% constraint: ignore certain components
idxcompts_ignore = mipopt.m_idxcompts_ignore;
tf = ismember(samplegraph.m_segid_per_compts',segids_forbid_always','rows');
idxcompts_forbid= find(tf);

lbsforbid = [GLOBALCONST.TABLE_MIDSUPPORT];
tf = ismember(samplegraph.m_complabels,lbsforbid);
idxcompts_forbid_by_label = find(tf);

idxcompts_forbid_by_label = [];

% also ignore components not covering specified seglbs
if ~isempty(mipopt.m_seglbs_use)
    idxcompts_ignore_by_seglbs = find(~any(seglbs2comp(mipopt.m_seglbs_use,:),1));
    
    % if we have fixed component, do not add them to ignore list
    if ~isempty(mipopt.m_idxcompts_fix)
        idxcompts_ignore_by_seglbs = setdiff(idxcompts_ignore_by_seglbs,mipopt.m_idxcompts_fix);
    end
else
    idxcompts_ignore_by_seglbs = [];
end
% idxcompts_ignore_by_seglbs = [];
idxcompts_ignore = cat(2,idxcompts_ignore,idxcompts_forbid(:)',idxcompts_ignore_by_seglbs(:)',idxcompts_forbid_by_label(:)');
idxcompts_ignore = unique(idxcompts_ignore(:))';

ijks = zeros(3,length(idxcompts_ignore));
if ~isempty(idxcompts_ignore)
    ijks(1,:) = 1:length(idxcompts_ignore);
    xx = x2varidx(idxcompts_ignore);
    ijks(2,:) = xx;
    ijks(3,:) = 1;
end
if ~isempty(ijks)
    coefAeq_ignorecompts = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
else
    coefAeq_ignorecompts = sparse(0,nvar);
end
beq_ignorecompts = zeros(size(coefAeq_ignorecompts,1),1);

%% hard-assignment definition: 
% each sketch segment is assigned to an element, and conflicting segments
% cannot assign to the same element
% v(i,j)==1 iff ith segment is assigned to jth element

% constraints and objectives:
% v(i,j)==0 iff ith segment is not covered by jth element (automatic guaranteed)
% v(i,j)<=x(j), v(i,j) can be hard-assigned to element j only if the component is selected
% v(k1,j)+v(k2,j)<=1 iff segment k1,k2 conflict at element j
% sum[i](v(i,j)*coverlen(i))>=ratio*x(j), the segments assigned to element j
% should cover large enough portion of the element, if the element is
% selected
ijks = zeros(3,0);
b_hardassign = [];
idxlastrow = 0;

ijks_eq = zeros(3,0);
idxlastrow_eq = 0;
beq_hardassign = [];

% coverlen benefit:
% encourage selecting elements of large coverlen
f_hardassign = zeros(nvar,1); 
for i=1:nx
    % seglbs covered by this component
    seglbs = find(seglbs2comp(:,i));
    xx = x2varidx(i);
    conflict_seglb2seglb = samplegraph.m_conflict_seg2seg2comp(:,:,i);
    [maxcoverlen,this_conflictgroup] = ...
        samplegraph.get_max_coverlen_per_comp(mipopt.m_covermode,i);
    cgroup = this_conflictgroup{1};
    
    % ignore conflicting segments that are not covered in current coverage
    maskvalid = false(size(conflict_seglb2seglb));
    maskvalid(seglbs,seglbs) = true;
    conflict_seglb2seglb(~maskvalid) = false;
    
    % filtering based on whether assignment is use
    if mipopt.m_usesegassign
        % dynamic assignment is used, so segment conflict is taken into
        % consideration
        
        % conflict: v1+v2<=1
        [s1,s2] = find(tril(conflict_seglb2seglb));
        v1 = idxmap_seglbs2comp(s1,i);
        v2 = idxmap_seglbs2comp(s2,i);
        vv1 = v2varidx(v1);
        vv2 = v2varidx(v2);
        ents = zeros(3,length(v1)*2);
        ents(1,:) = [1:length(v1),1:length(v1)] + idxlastrow;
        ents(2,:) = [vv1(:)',vv2(:)'];
        ents(3,:) = 1;
        ijks = cat(2,ijks,ents);
        idxlastrow = idxlastrow + length(v1);
        b_tmp = ones(length(v1),1);
        b_hardassign = cat(1,b_hardassign,b_tmp);
        
        % conflict group advantage prevention:
        % for each confliction group, one seglbs must be selected, preventing
        % favoring element with lots of conflicts because it requires less
        % segment assignement and generates less overlap
        % sum[k](v(k))>=x, iff v(k) in the same conflict group
        for k=1:length(cgroup)
            cseglbs = cgroup{k};
            tf = ismember(seglbs,cseglbs);
            if nnz(tf) <= 1
                %no conflict, ignore
                continue;
            end
            
            seglbs_incgroup = seglbs(tf);
            vs = idxmap_seglbs2comp(seglbs_incgroup,i);
            vv = v2varidx(vs);
            ents = zeros(3,length(vv)+1);
            ents(1,:) = idxlastrow + 1;
            ents(2,:) = [xx,vv(:)'];
            ents(3,:) = [1,-ones(1,length(vv))];
            b_hardassign(end+1) = 0;
            ijks = cat(2,ijks,ents);
            idxlastrow = idxlastrow + 1;
        end
    end
    
    % fix value: if a segment has no conflict, it must be hard-assigned to
    % this component if the component is selected
    % v(i,j) = x(j) iff v(i,j) has no conflict
    if mipopt.m_usesegassign
        % dynamic assignment, pre-assign only for conflict-free segments
        s3 = find(any(conflict_seglb2seglb,1));
        seglbs_noconflict = setdiff(seglbs,s3);
    else
        % static assignment, assign if the segment is covered
        seglbs_noconflict = seglbs;
    end
    vs = idxmap_seglbs2comp(seglbs_noconflict,i);
    vv = v2varidx(vs);
    fixvars(vv) = true;
    ents = zeros(6,length(vv));
    ents(1,:) = idxlastrow_eq + (1:length(vv));
    ents(2,:) = vv;
    ents(3,:) = 1;
    ents(4,:) = ents(1,:);
    ents(5,:) = xx;
    ents(6,:) = -1;
    ijks_eq = cat(2,ijks_eq,ents(1:3,:),ents(4:6,:));
    idxlastrow_eq = idxlastrow_eq + length(vv);
    
    % continuous group integrity:
    % if a bunch of segments belong to a contiuous group, then they are
    % simultaneously assigned or not assigned to a component
    this_contgroup2seglbs = false(size(contgroup2seglbs));
    this_contgroup2seglbs(:,seglbs_noconflict) = contgroup2seglbs(:,seglbs_noconflict); %ignore seglbs not involved in this component
    for k=1:size(this_contgroup2seglbs,1)
        seglbs_ingroup = find(this_contgroup2seglbs(k,:));
        vs = idxmap_seglbs2comp(seglbs_ingroup,i);
        vv = v2varidx(vs);
        if length(vv)>1
            % v(k)=v(1) for all v(k)
            nrow = length(vv)-1;
            ents = zeros(6,nrow);
            ents(1,:) = (1:nrow) + idxlastrow_eq;
            ents(2,:) = vv(1);
            ents(3,:) = -1;
            ents(4,:) = ents(1,:);
            ents(5,:) = vv(2:end);
            ents(6,:) = 1;
            ijks_eq = cat(2,ijks_eq,ents(1:3,:),ents(4:6,:));
            idxlastrow_eq = idxlastrow_eq + nrow;
        end
    end
    
    % sufficient coverage constraint:
    % sum[i](v(i,j)*coverlen(i))>=ratio*x(j)
%     vs = idxmap_seglbs2comp(seglbs,i);
%     vv = v2varidx(vs);
%     vlens = samplegraph.m_coverlen_seg2comp(seglbs,i);
% %     maxcoverlen = sum(samplegraph.m_coverlen_seg2comp(seglbs,i));
%     maxcoverlen = max(1,maxcoverlen) * 0.7;
%     ents = zeros(3,length(vs)+1);
%     ents(1,:) = idxlastrow + 1;
%     ents(2,:) = [xx, vv(:)'];
%     ents(3,:) = [maxcoverlen,-ones(1,length(vv)).*(vlens(:)')];
%     ijks = cat(2,ijks,ents);
%     idxlastrow = idxlastrow+1;
%     b_hardassign(end+1) = 0;
    
    % validity: v(i,j)<=x(j)
    vs = idxmap_seglbs2comp(seglbs,i);
    vv = v2varidx(vs);
    ents = zeros(6,length(vv));
    ents(1,:) = (1:length(vv))+idxlastrow;
    ents(2,:) = vv;
    ents(3,:) = 1;
    ents(4,:) = ents(1,:);
    ents(5,:) = xx;
    ents(6,:) = -1;
    ijks = cat(2,ijks,ents(1:3,:),ents(4:6,:));
    idxlastrow = idxlastrow + length(vv);
    b_tmp = zeros(1,length(vv));
    b_hardassign = cat(1,b_hardassign(:),b_tmp(:));
    
    % objective: maximize sum[i](v(i,j)*len(i))-a*x(j)
    % tend to select elements that have higher coverlen, but each element
    % introduces a coverlen threshold a, if the total coverlen does not
    % reach a, generates a penalty
    vs = idxmap_seglbs2comp(seglbs,i);
    vv = v2varidx(vs);
    vlens = samplegraph.m_coverlen_seg2comp(seglbs,i);
%     vfitness = samplegraph.m_coverfitness_seg2comp(seglbs,i);
    f_hardassign(vv) = -vlens;
    f_hardassign(xx) = thres_min_coverlen_to_award;
    f_hardassign(xx) = max(1,maxcoverlen);
%     f_hardassign(:) = 0;
    if ~mipopt.m_usesegassign
        % not using dynamic assignment
        f_hardassign(:) = 0;
    end
end
if ~isempty(ijks)
    coefA_hardassign = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
else
    coefA_hardassign = sparse(0,nvar);
end
b_hardassign = b_hardassign(:);
if ~isempty(ijks_eq)
    coefAeq_hardassign = sparse(ijks_eq(1,:),ijks_eq(2,:),ijks_eq(3,:),max(ijks_eq(1,:)),nvar);
else
    coefAeq_hardassign = sparse(0,nvar);
end
beq_hardassign = zeros(size(coefAeq_hardassign,1),1);

%% coverage benefit
u2covercost = -u2seglens/nsample;

%% match cost
compscores = mat2gray(samplegraph.m_compscores);
if nbin_match >0
    for i=1:length(samplegraph.m_dblbs)
        thislb = samplegraph.m_dblbs(i);
        mask = samplegraph.m_complabels == thislb;
        thiscores = compscores(mask);
        if isempty(thiscores)
            continue;
        end
        thiscores = mat2gray(thiscores);
        thiscores = util_quantize(thiscores,nbin_match);
        compscores(mask) = thiscores;
    end
end
matchcost_per_comp = 1-compscores;
x2matchcost = matchcost_per_comp;

%% structure cost
f_struct = zeros(nvar,1);
yy = y2varidx(1:ny);
f_struct(yy) = y2structcost;

%% overlap cost
f_overlap = zeros(nvar,1);
cc = c2varidx(1:nc);
f_overlap(cc) = c2ovcost;

%% overlap definition
% c(k) = sum(v(k,:))-u(k)
ijks = zeros(3,0);
beq_overlap = [];
for i=1:nc
    % ith 2d segment are covered by which components?
    vs = seglbs2vs(i);
    
    vv = v2varidx(vs);
    cc = c2varidx(i);
    uu = u2varidx(i);
    ents = zeros(3,2+length(vv));
    
    % ci+ui-sum(x)=0
    ents(1,:) = i;
    ents(2,:) = [cc,uu,vv(:)'];
    ents(3,:) = [1,1,-ones(1,length(vv))];
    ijks = cat(2,ijks,ents);
    
    if ismember(i,seglbs_always_covered)
        % ci+ui-sum(x)=1 --> ci=sum(x)
        beq_overlap(end+1) = 1;
    else
        beq_overlap(end+1) = 0;
    end
end
coefAeq_overlap = sparse(ijks(1,:),ijks(2,:),ijks(3,:),nc,nvar);
% beq_overlap = zeros(nc,1);
beq_overlap = beq_overlap(:);

%% coverage definition and objective
ijks = zeros(3,0);
idxlastrow = 0;
b_cover = [];
for i=1:nu
    % hard-assignment variables that refer to this segment
    vs = seglbs2vs(i);
    vv = v2varidx(vs);
    uu = u2varidx(i);
    
    if ismember(i,seglbs_always_covered)
        % this seglbs is covered by ghost component
        % u>=1 --> -u<=-1
        ents = [idxlastrow+1,uu,-1]';
        ijks = cat(2,ijks,ents);
        idxlastrow = idxlastrow+1;
        b_cover(end+1) = -1;
    else
        % u<=sum[k](v(k))
        ents = zeros(3,1+length(vs));
        ents(1,:) = idxlastrow+1;
        ents(2,:) = [uu,vv(:)'];
        ents(3,:) = [1,-ones(1,length(vv))];
        idxlastrow = idxlastrow+1;
        ijks = cat(2,ijks,ents);
        b_cover(end+1) = 0;

        % u>=v(k)
        ents = zeros(3,2*length(vs));
        for k=1:length(vs)
            t = idxlastrow+1;
            idxlastrow = idxlastrow+1;
            ents(:,2*k-1) = [t,vv(k),1];
            ents(:,2*k) = [t,uu,-1];
            b_cover(end+1) = 0;
        end
        ijks = cat(2,ijks,ents);
    end
end
coefA_cover = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
% b_cover = zeros(size(coefA_cover,1),1);
b_cover = b_cover(:);
f_cover = zeros(nvar,1);
uu = u2varidx(1:nu);
f_cover(uu) = u2covercost;

%% minimum cover constraint
uu = u2varidx(1:nu);
coefA_cover(end+1,uu) = -u2seglens;
if ~isempty(mipopt.m_mincover_in_seglbsuse)
    lens_always_cover = sum(u2seglens(seglbs_always_covered));
    lens_must_cover  = lens_always_cover + (sum(u2seglens)-lens_always_cover)*mipopt.m_mincover_in_seglbsuse;
else
    lens_must_cover = sum(u2seglens)*mipopt.m_mincover;
end
b_cover(end+1) = -lens_must_cover;
% b_cover(end+1) = -sum(u2seglens) * mipopt.m_mincover; %at least cover this much

%% match cost
f_match = zeros(nvar,1);
xx = x2varidx(1:nx);
f_match(xx) = x2matchcost;

%% symmetric constraint
x=1:nx;
xsym = samplegraph.get_symmetric_idxcompts(x);
maskhavesym = xsym>0;
x = x(maskhavesym);
xsym = xsym(maskhavesym);
sympairs = cat(1,x(:)',xsym(:)');
sympairs = sort(sympairs,1);
sympairs = unique(sympairs','rows')';
ijks = zeros(3,0);

if mipopt.m_usesym
    for i=1:size(sympairs,2)
        x1 = sympairs(1,i);
        x2 = sympairs(2,i);
        xx1 = x2varidx(x1);
        xx2 = x2varidx(x2);
        ijks(:,end+1) = [i,xx1,1];
        ijks(:,end+1) = [i,xx2,-1];
    end
end
if ~isempty(ijks)
    coefAeq_symmetry = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
else
    coefAeq_symmetry = sparse(0,nvar);
end
beq_symmetry = zeros(size(coefAeq_symmetry,1),1);

% degree constraint
dblbs = samplegraph.m_dblbs;
maxdgmat_lbs2lbs = zeros(max(dblbs));
maxdgmat_lbs2lbs(dblbs,dblbs) = samplegraph.m_maxdgmat;

if GLOBALCONST.DBG_CHECK('connect if cover','MIP: dropping max degree constraint')
    maxdgmat_lbs2lbs(:) = 10;
end

edgeidxmap = double(tril(adjmat));
edgeidxmap(edgeidxmap>0) = 1:nnz(edgeidxmap);
edgeidxmap = edgeidxmap + edgeidxmap';
b_maxdg = [];
idxlastrow = 0;
ijks = zeros(3,0);
for i=1:length(dblbs)
    lbfrom = dblbs(i);
    
    % components that belong to this label
    idxcompts_lbfrom = find(samplegraph.m_complabels==lbfrom);
    if isempty(idxcompts_lbfrom)
        continue;
    end
    
    for j=1:length(dblbs)
        lbto = dblbs(j);
        maxdg = maxdgmat_lbs2lbs(lbfrom,lbto); %lbfrom can have this many lbto
        idxcompts_lbto = find(samplegraph.m_complabels == lbto);
        if isempty(idxcompts_lbto)
            continue;
        end
        
        for k = 1:length(idxcompts_lbfrom)
            y = edgeidxmap(idxcompts_lbfrom(k),idxcompts_lbto);
            y = y(y~=0);
            yy = y2varidx(y);
            
            if length(y) <= maxdg
                % this constraint will not be violated anyway
                continue;
            end
            
            ents = zeros(3,length(yy));
            ents(1,:) = idxlastrow+1;
            ents(2,:) = yy;
            ents(3,:) = 1;
            ijks = cat(2,ijks,ents);
            b_maxdg(end+1) = maxdg;
            
            idxlastrow = idxlastrow+1;
        end
    end
end
if ~isempty(ijks)
    coefA_maxdg = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
else
    coefA_maxdg = sparse(0,nvar);
end
b_maxdg = b_maxdg(:);

%% label count constraint: number of components of a particular label must
% be within a range
lbs2maxcount = GLOBALCONST.LABEL_MAXCOUNT;
lbs2mincount = GLOBALCONST.LABEL_MINCOUNT;

if ~isempty(mipopt.m_labelcount)
    lbscount = mipopt.m_labelcount;
    if isvector(lbscount)
        lbscount = lbscount(:);
    end
    for i=1:size(lbscount,2)
        thislb = lbscount(1,i);
        mincount = lbscount(2,i);
        maxcount = lbscount(3,i);
        lbs2maxcount(thislb) = min(lbs2maxcount(thislb),maxcount);
        lbs2mincount(thislb) = max(lbs2mincount(thislb),mincount);
    end
end

ijks = zeros(3,0);
dblbs = samplegraph.m_dblbs;
b_labelcount = [];
idxlastrow = 0;
for i=1:length(dblbs)
    maxcount_thislb = lbs2maxcount(dblbs(i));
    mincount_thislb = lbs2mincount(dblbs(i));
    
    x = find(samplegraph.m_complabels == dblbs(i));
    if isempty(x)
        continue;
    end
    
    xx = x2varidx(x);
    
    if maxcount_thislb~=inf
        % max count: sum(xk)<=m
        ents = zeros(3,length(x));
        ents(1,:) = idxlastrow+1;
        ents(2,:) = xx;
        ents(3,:) = 1;
        ijks = cat(2,ijks,ents);
        idxlastrow = idxlastrow+1;
        b_labelcount(end+1) = maxcount_thislb;
    end
    
    if mincount_thislb >0
        % min count: sum(xk)>=e
        ents = zeros(3,length(x));
        ents(1,:) = idxlastrow+1;
        ents(2,:) = xx;
        ents(3,:) = -1;
        ijks = cat(2,ijks,ents);
        idxlastrow = idxlastrow+1;
        b_labelcount(end+1) = -mincount_thislb;
    end
end
if ~isempty(ijks)
    coefA_labelcount = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
else
    coefA_labelcount = sparse(0,nvar);
end
b_labelcount = b_labelcount(:);

%% total component count constraint
mincomp = mipopt.m_mincomp;
maxcomp = mipopt.m_maxcomp;
ijks = zeros(3,0);
b_compcount = [];
idxlastrow = 0;
x = 1:nx;
xx = x2varidx(x);
if mincomp > 0
    ents = zeros(3,length(xx));
    ents(1,:) = idxlastrow+1;
    ents(2,:) = xx;
    ents(3,:) = -1;
    ijks = cat(2,ijks,ents);
    idxlastrow = idxlastrow+1;
    b_compcount(end+1) = -mincomp;
end

if maxcomp ~= inf
    ents = zeros(3,length(xx));
    ents(1,:) = idxlastrow+1;
    ents(2,:) = xx;
    ents(3,:) = 1;
    ijks = cat(2,ijks,ents);
    idxlastrow = idxlastrow+1;
    b_compcount(end+1) = maxcomp;
end
if ~isempty(ijks)
    coefA_compcount = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
else
    coefA_compcount = sparse(0,nvar);
end
b_compcount = b_compcount(:);

% fixed component constraint: some components are given
idxcompts_fix = mipopt.m_idxcompts_fix;
ijks = zeros(3,length(idxcompts_fix));
if ~isempty(idxcompts_fix)
    ijks(1,:) = 1:length(idxcompts_fix);
    ijks(2,:) = x2varidx(idxcompts_fix);
    ijks(3,:) = 1;
end
beq_fixcompts = ones(length(idxcompts_fix),1);
if ~isempty(ijks)
    coefAeq_fixcompts = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
else
    coefAeq_fixcompts = sparse(0,nvar);
end

%% unexposed/overlap definition (r variables) and constraint
if ~use_overlap_hard_cons
    % disable overlap hard constraint as requested
    coefA_unexpose = sparse(0,nvar);
    b_unexpose = [];
else
    ijks = zeros(3,0);
    b_unexpose = [];
    idxlastrow = 0;
    ncomp_per_seg = sum(seglbs2comp,2); %upper bound of c variable
    for i=1:nx
        x = i;
        rs = x2rs(x);
        us = x2us(x);
%         vs = seglbs2vs(us);
        vs = idxmap_seglbs2comp(us,i);
        
        if isempty(us)
            % a component does not cover any segment, ignore it
            continue;
        end

        rr = r2varidx(rs);
        cc = c2varidx(us); %cs == us, c and u are one to one corresponded
        vv = v2varidx(vs);

        % r(i)<=v(i), i.e. if segment u(i) is not hard-assigned to this
        % element, it can't be overlapped by others
        ents = zeros(6,length(rs)); %[i1,j1,k1 i2,j2,k2]
        ents(1,:) = idxlastrow+(1:length(rs));
        ents(2,:) = rr;
        ents(3,:) = 1;
        ents(4,:) = ents(1,:);
        ents(5,:) = vv;
        ents(6,:) = -1;
        ijks = cat(2,ijks,ents(1:3,:),ents(4:6,:));
        idxlastrow = idxlastrow+length(rs);
        b_unexpose(end+(1:length(rs))) = 0;

        % r(i)<=c(i), i.e. if segment u(i) is not hard-assigned twice, there's
        % not way for overlapping to happen
        ents = zeros(6,length(rs));
        ents(1,:) = idxlastrow+(1:length(rs));
        ents(2,:) = rr;
        ents(3,:) = 1;
        ents(4,:) = ents(1,:);
        ents(5,:) = cc;
        ents(6,:) = -1;
        ijks = cat(2,ijks,ents(1:3,:),ents(4:6,:));
        idxlastrow = idxlastrow+length(rs);
        b_unexpose(end+(1:length(rs))) = 0;

        % r(i)>=v(i)+c/cmax-1, only if segment u(i) is hard-assigned to
        % this element, and it's also assigned to others, can overlapping
        % hapen
        cmax = ncomp_per_seg(us); %max possible c
        cmax(cmax==0) = 1;
        ents = zeros(9,length(rs));
        ents(1,:) = idxlastrow+(1:length(rs));
        ents(2,:) = vv;
        ents(3,:) = 1;
        ents(4,:) = ents(1,:);
        ents(5,:) = rr;
        ents(6,:) = -1;
        ents(7,:) = ents(1,:);
        ents(8,:) = cc;
        ents(9,:) = 1./cmax;
        ijks = cat(2,ijks,ents(1:3,:),ents(4:6,:),ents(7:9,:));
        idxlastrow = idxlastrow+length(rs);
        b_unexpose(end+(1:length(rs))) = 1;

        % sum(rs.*lens) <= max_overlap_ratio * sum(lens)
        rlens = u2seglens(us);
        ents = zeros(3,length(rs));
        ents(1,:) = idxlastrow+1;
        ents(2,:) = rr;
        ents(3,:) = rlens;
        ijks = cat(2,ijks,ents);
        idxlastrow = idxlastrow+1;
        b_unexpose(end+1) = max_overlap_ratio * sum(rlens);
    end
    coefA_unexpose = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
    b_unexpose = b_unexpose(:);
end

%% 2d segment cover constraint: long segments must be covered
% clbs = clusterdata(u2seglens(:),3);
kcenter = [max(u2seglens),mean(u2seglens),median(u2seglens),min(u2seglens)]';
clbs = kmeans(u2seglens(:),4,'emptyaction','singleton','start',kcenter);
lens_per_cluster = zeros(1,4);
for i=1:4
    lens_per_cluster(i) = mean(u2seglens(clbs==i));
end
[~,idxsort] = sort(lens_per_cluster,'descend');
idxuse = idxsort(1:2);
maskuse = ismember(clbs,idxuse);
% maskuse = u2seglens > mean(u2seglens);
maskcoverable = any(seglbs2comp,2);
maskuse = maskuse(:)' & maskcoverable(:)';
us = find(maskuse);
uu = u2varidx(us);
f_mustcover = zeros(nvar,1);
f_mustcover(uu) = -u2seglens(us)/sum(u2seglens);
% ijks = zeros(3,length(us));
% ijks(1,:) = 1:length(us);
% ijks(2,:) = uu;
% ijks(3,:) = 1;
% beq_mustcover = ones(length(us),1);
% if ~isempty(ijks)
%     coefAeq_mustcover = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
% else
%     coefAeq_mustcover = sparse(0,nvar);
% end

%% define variables h: h(i,j)=1 iff ith segment is assigned jth label
% h(i,j)>=v(k), v(k) is segment i hard-assigned to an element of label j
% h(i,j)<=sum(v(k)), combined with above, equivalent to 'or'
ijks = zeros(3,0);
idxlastrow = 0;
for i=1:nh
    h = i;
    seglbs = h2seglabelpair(1,i);
    label = h2seglabelpair(2,i);
    vs = seglbs_label2vs(seglbs,label);
    
    % h>=v(k)
    hh = h2varidx(h);
    vv = v2varidx(vs);
    ents = zeros(6,length(vv));
    ents(1,:) = idxlastrow + (1:length(vv));
    idxlastrow = idxlastrow + length(vv);
    ents(2,:) = vv;
    ents(3,:) = 1;
    ents(4,:) = ents(1,:);
    ents(5,:) = hh;
    ents(6,:) = -1;
    ijks = cat(2,ijks,ents(1:3,:),ents(4:end,:));
    
    % h<=sum(v(k))
    ents = zeros(3,length(vv)+1);
    ents(1,:) = idxlastrow + 1;
    idxlastrow = idxlastrow+1;
    ents(2,:) = [hh,vv(:)'];
    ents(3,:) = [1,-ones(1,length(vv))];
    ijks = cat(2,ijks,ents);
end
if ~isempty(ijks)
    coefA_seglabel = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
else
    coefA_seglabel = sparse(0,nvar);
end
b_seglabel = zeros(size(coefA_seglabel,1),1);

%% segment group constraint/cost
% some 2d segments preferrably take the same label
% t(i,j)>=h(k,j) for all segment k that belongs to group i that end up taking label j
% t(i,j)<=sum(h(k,j)), if no segment in this group take label j, t=0
ijks = zeros(3,0);
idxlastrow = 0;
f_segroup = zeros(nvar,1);
totalseglens = sum(lens_per_2dseg);
for i=1:nt
    % group id for this t
    gid = t2grouplabelpair(1,i);
    
    % intended label for this group
    label = t2grouplabelpair(2,i);
    
    % segments in this group
    seglbs = seglbsgroup{gid};
    seglens = lens_per_2dseg(seglbs);
    
    % find h(i) that relates to (seglbs(i),label(i))
    lbs = repmat(label,1,length(seglbs));
    idx = sub2ind(size(hidx_seglbs2label),seglbs(:),lbs(:));
    hs = hidx_seglbs2label(idx);
    
    % sometimes you want segment (a,b,c,d) to be covered by a particularly 
    % type of component but the matching does not find any of these components
    % covering one of the segments, rule out these requirements
    hs = hs(hs>0); %rule out impossible h
    seglens_per_h = seglens(hs>0);
    if isempty(hs)
        % this group will never take the label, ignore it
        continue;
    end
    
    % t>=h(k)
    tt = t2varidx(i);
    hhs = h2varidx(hs);
    ents = zeros(6,length(hhs));
    ents(1,:) = (1:length(hhs))+idxlastrow;
    idxlastrow = idxlastrow+length(hhs);
    ents(2,:) = hhs;
    ents(3,:) = 1;
    ents(4,:) = ents(1,:);
    ents(5,:) = tt;
    ents(6,:) = -1;
    ijks = cat(2,ijks,ents(1:3,:),ents(4:end,:));
    
    % t<=sum(h(k))
    ents = zeros(3,length(hhs)+1);
    ents(1,:) = idxlastrow + 1;
    idxlastrow = idxlastrow+1;
    ents(2,:) = [tt,hhs(:)'];
    ents(3,:) = [1,-ones(1,length(hhs))];
    ijks = cat(2,ijks,ents);
    
    % cost function
    % f=kt-sum(h(i)), k=number of segments in this group
    % when all(h(i)==0), t=0, f=0
    % when any(h(i)==1), t=1, f=k-sum(h(i)), only vanish when all h=1
    k = length(seglbs);
    % f_segroup(tt) = f_segroup(tt) + k * sum(seglens)/totalseglens;
    
    % in the above formulation, each adverse segment induces cost 1, here
    % we refine it, each adverse segment induces cost |seglength|, so the
    % "total length" k becomes the real total length sum(seglens).
    f_segroup(tt) = f_segroup(tt) + sum(seglens)/totalseglens;
    f_segroup(hhs) = f_segroup(hhs) - 1 * seglens_per_h(:)/totalseglens;
end
if ~isempty(ijks)
    coefA_segroup = sparse(ijks(1,:),ijks(2,:),ijks(3,:),max(ijks(1,:)),nvar);
else
    coefA_segroup = sparse(0,nvar);
end
b_segroup = zeros(size(coefA_segroup,1),1);

%% combine model
coefA_opt = cat(1,coefA_cover,coefA_labelcount,...
    coefA_maxdg,coefA_compcount,coefA_unexpose,...
    coefA_seglabel,coefA_segroup,coefA_hardassign);

b_opt = cat(1,b_cover,b_labelcount,...
    b_maxdg,b_compcount,b_unexpose,...
    b_seglabel,b_segroup,b_hardassign);

coefAeq_opt = cat(1,coefAeq_overlap,coefAeq_symmetry,...
    coefAeq_fixcompts,coefAeq_hardassign,coefAeq_ignorecompts);
beq_opt = cat(1,beq_overlap,beq_symmetry,...
    beq_fixcompts,beq_hardassign,beq_ignorecompts);
% f_opt = cat(1,f_cover(:),f_match(:),f_overlap(:),f_struct(:));

f_opt = f_cover * w_cover + ...
    f_mustcover*w_mustcover + ...
    f_match*w_match + ...
    f_overlap*w_overlap + ...
    f_struct*w_struct + ...
    f_segroup*w_segroup + ...
    f_hardassign*w_hardassign;

% add graph connectivity constraint
[gA,gb,gAeq,gbeq,gidxbin,sol2part,sos] = util_formulate_graph_connect(adjmat | adjmat');

% combine the 2
nhead = nx+ny;
nmid = size(gA,2) - nx - ny;
ntail = nvar - nx-ny;
nvar = nhead + nmid + ntail; %new variable count after adding flow-related variables
coefA = sparse(size(coefA_opt,1),nvar);
coefA(:,1:nhead) = coefA_opt(:,1:nhead);
coefA(:,end-ntail+1:end) = coefA_opt(:,end-ntail+1:end);
b = b_opt;
coefAeq = sparse(size(coefAeq_opt,1),nvar);
coefAeq(:,1:nhead) = coefAeq_opt(:,1:nhead);
coefAeq(:,end-ntail+1:end) = coefAeq_opt(:,end-ntail+1:end);
beq = beq_opt;
f = zeros(nvar,1);
f(1:nhead) = f_opt(1:nhead);
f(end-ntail+1:end) = f_opt(end-ntail+1:end);

gA(:,end+ntail) = 0;
gAeq(:,end+ntail) = 0;
coefA = cat(1,coefA,gA);
b = cat(1,b,gb);
coefAeq = cat(1,coefAeq,gAeq);
beq = cat(1,beq,gbeq);

% be careful here, re-index variables
idxbin = cat(1,gidxbin(:),...
    nhead+nmid+(1:nu)',...
    nhead+nmid+nu+nc+(1:nr)',...
    nhead+nmid+nu+nc+nr+(1:nh)',...
    nhead+nmid+nu+nc+nr+nh+(1:nt)',...
    nhead+nmid+nu+nc+nr+nh+nt+(1:nv)');

fixvars = fixvars(:);
fixvars = cat(1,fixvars(1:nhead),nan(nmid,1),fixvars(end-ntail+1:end));

% delete fixed variables in the model
% idxvarfixed = find(~isnan(fixvars));
% idxvarremain = find(isnan(fixvars));
% valuefixed = fixvars(idxvarfixed);
% f_trim = f(idxvarremain);
% 
% coefA_fixed = coefA(:,idxvarfixed);
% coefA_trim = coefA(:,idxvarremain);
% b_trim = b - coefA_fixed*valuefixed(:);
% 
% coefAeq_fixed = coefAeq(:,idxvarfixed);
% coefAeq_trim = coefAeq(:,idxvarremain);
% beq_trim = beq - coefAeq_fixed*valuefixed(:);
% 
% maskbin = false(nvar,1);
% maskbin(idxbin) = true;
% maskbin = maskbin(idxvarremain);
% idxbin_trim = find(maskbin);
% n_fixed_before_sos = nnz(idxvarfixed < min(sos.index));
% sos.index = sos.index - n_fixed_before_sos;
% gurobi_model = util_gurobi_from_mip(f_trim,coefA_trim,b_trim,coefAeq_trim,beq_trim,idxbin_trim,[],sos);

gurobi_model = util_gurobi_from_mip(f,coefA,b,coefAeq,beq,idxbin,[],sos);
linmodel.fixvars = fixvars;
linmodel.nvar_at_formulation = size(coefA,2);
linmodel.A = coefA;
linmodel.b = b;
linmodel.Aeq = coefAeq;
linmodel.beq = beq;
linmodel.v2varidx = v2varidx + nmid;
linmodel.y2varidx = y2varidx;
linmodel.f = f;

linmodel.f_cover = zeros(nvar,1);
linmodel.f_cover([1:nhead,end-ntail+1:end]) = f_cover*w_cover;

linmodel.f_mustcover = zeros(nvar,1);
linmodel.f_mustcover([1:nhead,end-ntail+1:end]) = f_mustcover*w_mustcover;

linmodel.f_match = zeros(nvar,1);
linmodel.f_match([1:nhead,end-ntail+1:end]) = f_match*w_match;

linmodel.f_overlap = zeros(nvar,1);
linmodel.f_overlap([1:nhead,end-ntail+1:end]) = f_overlap*w_overlap;

linmodel.f_struct = zeros(nvar,1);
linmodel.f_struct([1:nhead,end-ntail+1:end]) = f_struct*w_struct;

linmodel.f_segroup = zeros(nvar,1);
linmodel.f_segroup([1:nhead,end-ntail+1:end]) = f_segroup*w_segroup;

linmodel.f_hardassign = zeros(nvar,1);
linmodel.f_hardassign([1:nhead,end-ntail+1:end]) = f_hardassign*w_hardassign;

linmodel.idxbin = idxbin;
linmodel.sol2part = sol2part;
linmodel.nx = nx;
linmodel.ny = ny;
linmodel.nf = 2*ny;
linmodel.ng = nx;
linmodel.nu = nu;
linmodel.nc = nc;
linmodel.nr = nr;
linmodel.nh = nh;
linmodel.nt = nt;
linmodel.nv = nv;
linmodel.adjmat = adjmat | adjmat';
linmodel.ridx_seglbs2comp = idxmap_seglbs2comp;
linmodel.hidx_seglbs2label = hidx_seglbs2label;
linmodel.tidx_groupid2label = tidx_groupid2label;
linmodel.explain_solution = @(x,linmodel)explain_solution(x,linmodel);
linmodel.termlist = {}; %additional terms added later
linmodel.seglbs = lens_per_2dseg;

edgeidxmap = zeros(size(adjmat));
edgeidxmap(tril(adjmat)) = 1:nnz(tril(adjmat));
edgeidxmap = edgeidxmap + edgeidxmap';
linmodel.edgeidxmap = edgeidxmap; %used for y variable indexing

gurobi_model = srec_addcons_special(gurobi_model,linmodel,samplegraph);

if ~isempty(mipopt.m_idxcomptslist_per_combination_forbid)
    % forbid certain combinations
    idxcomptstlistlist = mipopt.m_idxcomptslist_per_combination_forbid;
    for i=1:length(idxcomptstlistlist)
        idxcomptslist = idxcomptstlistlist{i};
        gurobi_model = srec_addcons_forbid_combination(idxcomptslist,gurobi_model,linmodel);
    end
end

% % additional constraints
% if ~isempty(mipopt.m_idxcompts_per_cluster)
%    [gurobi_model,linmodel] = srec_addcons_clustering(gurobi_model,linmodel,mipopt.m_idxcompts_per_cluster);
% end

end

function solpart = explain_solution(linx,linmodel)
% extract detail solution from MIP solution

%the model may have additional variables due to constraints added after
%formulation
linx_original = linx;
linx = linx(1:linmodel.nvar_at_formulation); 
thres = 1e-6;
solpart = linmodel;

% fixvars = linmodel.fixvars;
% linxnew = fixvars;
% linxnew(isnan(fixvars)) = linx;
% linx = linxnew;

nhead = linmodel.nx+linmodel.ny+linmodel.nf+linmodel.ng;
sol_graph = linmodel.sol2part(linx(1:nhead));
namelist = fieldnames(sol_graph);
for i=1:length(namelist)
    solpart.(namelist{i}) = sol_graph.(namelist{i});
end
solpart.linx = linx_original;

% segment coverage: whether a segment is covered or not
nu = linmodel.nu;
u = linx(nhead+(1:nu));
masksegcover = abs(u)>thres;
solpart.masksegcover = masksegcover;

% segment cover count: how many times is segment i covered?
nc = linmodel.nc;
c = linx(nhead+nu+(1:nc));
c = round(c) + double(masksegcover);
solpart.segcovercount = c;

% segment overlap indicator: is the ith segment of jth object overlapped by
% others?
nr = linmodel.nr;
if nr>0
    r = linx(nhead+nu+nc+(1:nr));
    ridx_seglbs2comp = linmodel.ridx_seglbs2comp;
    maskoverlap_seglbs2comp = false(size(ridx_seglbs2comp));
    inds = find(ridx_seglbs2comp);
    ridx = ridx_seglbs2comp(inds);
    [ridx,idxsort] = sort(ridx);
    inds = inds(idxsort);
    maskoverlap_seglbs2comp(inds) = abs(r)>thres;
    solpart.maskoverlap_seglbs2comp = maskoverlap_seglbs2comp;
else
    solpart.maskoverlap_seglbs2comp = [];
end

% label assignment mask, does ith segment take jth label?
nh = linmodel.nh;
if nh>0
    h = linx(nhead+nu+nc+nr+(1:nh));
    hidx_seglbs2label = linmodel.hidx_seglbs2label;
    maskcover_seglbs2label = false(size(hidx_seglbs2label));
    inds = find(hidx_seglbs2label);
    hidx = hidx_seglbs2label(inds);
    [hidx,idxsort] = sort(hidx);
    inds = inds(idxsort);
    maskcover_seglbs2label(inds) = abs(h)>thres;
    solpart.maskcover_seglbs2label = maskcover_seglbs2label;
else
    solpart.maskcover_seglbs2label = [];
end

% group label activation
% mask(i,j)=1 iff some segment in group i is covered by label j
nt = linmodel.nt;
if nt>0
    t = linx(nhead+nu+nc+nr+nh+(1:nt));
    tidx_groupid2label = linmodel.tidx_groupid2label;
    maskcover_groupid2label = false(size(tidx_groupid2label));
    inds = find(tidx_groupid2label);
    tidx = tidx_groupid2label(inds);
    [tidx,idxsort] = sort(tidx);
    inds = inds(idxsort);
    maskcover_groupid2label(inds) = abs(t)>thres;
    solpart.maskcover_groupid2label = maskcover_groupid2label;
else
    solpart.maskcover_groupid2label = [];
end

% 2d segment hard assignment
nv = linmodel.nv;
idxmap_seglbs2comp = linmodel.ridx_seglbs2comp;
if nv>0
    v = linx(nhead+nu+nc+nr+nh+nt+(1:nv));
    idx = find(idxmap_seglbs2comp);
    values = idxmap_seglbs2comp(idx);
    [~,idxsort] = sort(values);
    idx = idx(idxsort);
    maskassign_seglbs2comp = false(size(idxmap_seglbs2comp));
    maskassign_seglbs2comp(idx) = abs(v)>thres;
    solpart.maskassign_seglbs2comp = maskassign_seglbs2comp;
else
    solpart.maskassign_seglbs2comp = [];
end

namelist = fieldnames(linmodel);
for i=1:length(namelist)
    thisname = namelist{i};
    if length(thisname)>=2 && isequal(thisname(1:2),'f_')
        % this is a cost component
        costname = thisname(3:end);
        costname = ['cost_',costname];
        thiscost = dot(linx,linmodel.(thisname));
        solpart.(costname) = thiscost;
    end
end
solpart.cost_total = dot(linx,linmodel.f);

% solution that does not contain the flow-related variables
nhead = linmodel.nx+linmodel.ny;
ntail = nu+nc+nr+nh+nt+nv;
solpart.xnoflow = linx([1:nhead,end-ntail+1:end]);

for i=1:length(solpart.termlist)
    termdata = solpart.termlist{i};
    termdata = termdata.explain_term(solpart,termdata);
    solpart.termlist{i} = termdata;
end

end

function rs = x2rs_(x,idxmap_seglbs2comp)
    submat_ = idxmap_seglbs2comp(:,x);
    rs = submat_(submat_>0);
    rs = rs(:)';
end

function us = x2us_(x,seglbs2comp)
    [us,~] = find(seglbs2comp(:,x));
    us = us(:)';
end



%{
pp = inputParser;
pp.addParamValue('mincomp',0); %minimum number of components in output
pp.addParamValue('maxcomp',inf); %maximum number of components in output
pp.addParamValue('mincover',0.7); %minimum percentage of cover
pp.addParamValue('covermode','normal',@(x)any(strcmpi(x,{'normal','fullcover','justenough'})));
pp.addParamValue('useoverlap',true); %overlap hard constraint
pp.addParamValue('idxcompts_ignore',[]);%do not choose thses components

% seglbsgroup{i} = [seglb1,seglb2...seglbn], a set of seglbs that prefers
% to take the same label
% seglbsgroup can be 2*n cell, where seglbsgroup(:,i) = {seglbs,lb},
% indicating the specified segments prefer to take the label lb at the same
% time
pp.addParamValue('seglbsgroup',[]);

% segment continuous group
% g{i} = [seglb1,seglb2...], means seglb1,seglb2... belongs to a continuous
% stroke
pp.addParamValue('segcontgroup',[]);

% labelcount is 3*n matrix
% labelcount(:,i) = [label,mincount,maxcount], the component count of this
% label is between mincount and maxcount
pp.addParamValue('labelcount',[],@(x)isempty(x)||numel(x)==3||size(x,1)==3);

% given components
pp.addParamValue('idxcompts_fix',[]);
pp.parse(varargin{:});
%}