基于MeSC与交感神经作用关系的压力水平与白发模拟系统和压力规划系统
一个本人的Matlab项目，可用于根据压力水平模拟白发水平，并根据工作情况给出白发量最少的合理的压力规划。

% 细胞仿真
clc;
clear;
% Raw data
sl0 = [0 1 2 3 4];                                        % stress_level
MeSC = [7.5 10 5 2.5 1];                                  %No. of MeScs within per HF
whr1 = 0.4;                                               % max white hair rate 已知MeSC为1时白发率为40%
whr0 = 0;                                                 % whr when no stress 设无压力时白发率为0
Fr = 0.33;                                                % 女性的白发率
Mr = 0.29;                                                % 男性的白发率

% spline插值获得白发率
x = [MeSC(end) MeSC(1)];
y = [whr1 whr0];
whr = spline(x,y,MeSC);                                   % 利用已知的白发率与MeSC的关系插值

% Data visualization
figure()
y = [sl0;MeSC];
X = categorical({'Group1','Group2','Group3','Group4','Group5'});
b = bar(X,y);hold on
xtips1 = b(1).XEndPoints;
ytips1 = b(1).YEndPoints;
labels1 = string(b(1).YData);
text(xtips1,ytips1,labels1,'HorizontalAlignment','center',...
    'VerticalAlignment','bottom')
xtips2 = b(2).XEndPoints;
ytips2 = b(2).YEndPoints;
labels2 = string(b(2).YData);
text(xtips2,ytips2,labels2,'HorizontalAlignment','center',...
    'VerticalAlignment','bottom')
ylim([0,12]);
grid on;
title('Experiment data');
xlabel('test group');
colororder({'k','k'})
yyaxis left
ylabel('Stress level & MeSC level');
yyaxis right
ylabel('White hair rate');
plot(X,whr,'linewidth',2,'color','b');hold on;       % 利用压力水平与MeSC的对应关系，得到压力水平和白发率关系
scatter(X,whr);hold on;
ylim([-0.2,0.5]);
legend('stress level','MeSC level','White hair rate');

% Parameter setting
upr = 0.0001;                                             % 头发更新率 ???
nh = 100000;                                              % 头发数量，一般10万
nh_l = sqrt(nh);                                          % 生成图片的长度
Initial = 20;
%{
% Stress data input & simulation
% sli = input('每周的压力水平？(0-4)\n');
sli = input('每周的睡眠分数？(0-4)\n');
sli = ((100-sli)./100).*4;
sli = sli';
n = size(sli);
n0 = 1:1:n;                                               % 已知点
n1 = 1:n/(7*n(1)):n;                                      % 插值点
sl = spline(n0,sli,n1);                                   % 对压力水平插值
nsl = size(sl);
% Stress visualization
figure();
grid on;
yyaxis left
p1 = plot(1:nsl(2),sl,'linewidth',1.25,'color','b');      % 每日压力情况曲线
hold on;
stem(1:7:nsl(2),sli,'linewidth',1.25,'LineStyle','-.',...
     'MarkerFaceColor','red',...
     'MarkerEdgeColor','green');hold on;
title('Personal data');
xlabel('days');
ylabel('stress level');

% White hair rate simulation & visualization
whr2 = spline(sl0,whr,sl);                                % 根据输入的压力情况得到实际的每日白发率曲线
nwhr2 = size(whr2);
yyaxis right
p2 = plot(1:nsl(2),whr2,'linewidth',1.5,'color','r');hold on;
ylabel('white hair rate');
colororder({'b','r'})
legend([p1,p2],{'stress level','White hair rate','MeSC level'});

MeSC2 = spline(whr,MeSC,whr2);                            % 插值估计每日MeSC数量
% p3 = plot(1:nsl(2),MeSC2,'linewidth',1.5,'color','k');hold on;

%{
% 椭球头部模型
figure()
[x, y, z] = ellipsoid(0,0,0,3,4,4.5,fix(nh_l));
surf(x, y, -z)
c = copper;
c = flipud(c);
colormap(c);
axis equal
%}

% Hair simulation & visualization
hair = zeros(fix(nh_l),fix(nh_l))+Initial;                % fix取整,设置头发初始灰度
hair = uint8(hair);
wi = 1;                                                   % white i,变白了的头发的索引
bi = 1;                                                   % black i,变黑了的头发的索引
%F = zeros(1,nsl(2)*fix(upr*nh));
F = moviein(nsl(2)*fix(upr*nh));
j = 1;
wx = randi(fix(nh_l),1,fix(nh_l));
wy = randi(fix(nh_l),1,fix(nh_l));
for t = 1:1:nsl(2)                                        % 对每一天
    x = fix((nh_l-1)*rand(upr*nh,1)+1);                   % 求该天被更新的头发坐标集合，rand生成数组
    y = fix((nh_l-1)*rand(upr*nh,1)+1);
    for i = 1:fix(upr*nh)
        i_whr = whr2(t);
        if i_whr > 0
            hair(x(i),y(i)) = hair(x(i),y(i))+255*i_whr;
            wx(wi) = x(i);                                % 记录下变白了的头发
            wy(wi) = y(i);
            wi = wi+1;
        else
            if (hair(wx(bi),wy(bi))+255*i_whr)>=Initial
                hair(wx(bi),wy(bi)) = ...
                    hair(wx(bi),wy(bi))+255*i_whr;        % 此前变白的头发将优先恢复
                bi = bi+1;
            end
        end
        hair2 = mat2gray(hair);
        fig = figure();
        fig.Visible = 'off';
        subplot(121)
        imshow(hair2);
        subplot(122)
        Mf = imread('MeSC.png');
        imshow(Mf);
        hold on;
        N = MeSC2(t);
        ym = 829*(1-0.45).*rand(round(N),1);
        xm = ones(1,round(N)).*0.7*318;
        xmr= -20+(20+20)*rand(1,round(N));
        xm = xm+xmr;
        scatter(xm,ym,60,'g','filled');                   % 模拟MeSC的数量与分布
        
        F(j) = getframe(fig);
        j = j+1;
    end
end
figure()
%{
set(gcf,'unit','normalized',...                           % 本来打算设置figure大小
    'position',[0.1,0.1,0.56,0.42])
xlim([0,0.5]);
ylim([0,0.8]);
%}
movie(F,5);
%}

% Optimization
% 设对每一时刻，一个迫近（一周内）ddl的压力值为0.7，第二周的ddl压力值为0.3
% 对未来两周

% ddl input & stress level count
ddl1 = input('第一周ddl数目：');
ddl2 = input('第二周ddl数目：');
ddl3 = input('第三周ddl数目：');

sl1 = 0.7*ddl1+0.3*ddl2;
sl2 = 0.7*ddl2+0.3*ddl3;

% 多项式拟合求出sl-whr函数表达式
p=polyfit(sl0,whr,4);
whr1 = polyval(p,sl0);
figure()
plot(sl0,whr1,'linewidth',5,'color','b');hold on
plot(sl0,whr,'linewidth',2.5,'color','g');grid on
title('Polynomial fitting');
legend('Polynomial output','truth');
xlabel('stress level');
ylabel('white hair rate');


fun = @(x)-(7*p(1)+p(2)*(x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)+x(8)+x(9)+x(10)+x(11)+x(12)+x(13)+x(14))...
    +p(2)*(x(1)^2+x(2)^2+x(3)^2+x(4)^2+x(5)^2+x(6)^2+x(7)^2+x(8)^2+x(9)^2+x(10)^2+x(11)^2+x(12)^2+x(13)^2+x(14)^2)...
    +p(3)*(x(1)^3+x(2)^3+x(3)^3+x(4)^3+x(5)^3+x(6)^3+x(7)^3+x(8)^3+x(9)^3+x(10)^3+x(11)^3+x(12)^3+x(13)^3+x(14)^3)...
    +p(4)*(x(1)^4+x(2)^4+x(3)^4+x(4)^4+x(5)^4+x(6)^4+x(7)^4+x(8)^4+x(9)^4+x(10)^4+x(11)^4+x(12)^4+x(13)^4+x(14)^4));
x0 = [0,0,0,0,0,0,0,0,0,0,0,0,0,0];
A = [-1,-1,-1,-1,-1,-1,-1,0,0,0,0,0,0,0;...                 % 保证第一周ddl完成
    0,0,0,0,0,0,0,-1,-1,-1,-1,-1,-1,-1;...                  % 保证第二周ddl完成
    1,1,0,0,0,0,0,0,0,0,0,0,0,0;...                         % 别太累了，连续两天压力值不超过6
    0,1,1,0,0,0,0,0,0,0,0,0,0,0;...
    0,0,1,1,0,0,0,0,0,0,0,0,0,0;...
    0,0,0,1,1,0,0,0,0,0,0,0,0,0;...
    0,0,0,0,1,1,0,0,0,0,0,0,0,0;...
    0,0,0,0,0,1,1,0,0,0,0,0,0,0;...
    0,0,0,0,0,0,1,1,0,0,0,0,0,0;...
    0,0,0,0,0,0,0,1,1,0,0,0,0,0;...
    0,0,0,0,0,0,0,0,1,1,0,0,0,0;...
    0,0,0,0,0,0,0,0,0,1,1,0,0,0;...
    0,0,0,0,0,0,0,0,0,0,1,1,0,0;...
    0,0,0,0,0,0,0,0,0,0,0,1,1,0;...
    0,0,0,0,0,0,0,0,0,0,0,0,1,1;];                         % 不等式约束的系数矩阵，>号，要取相反数
b = [-sl1*7;-sl2*7;6;6;6;6;6;6;6;6;6;6;6;6;6];             % 不等式约束的常向量，>号，要取相反数
Aeq = [];                                                  % 等式约束的系数矩阵
beq = [];                                                  % 等式约束的常向量
lb = [0;0;0;0;0;0;0;0;0;0;0;0;0;0;];                       % 自变量的最小值
ub = [4;4;4;4;4;4;4;4;4;4;4;4;4;4;];                       % 自变量的最大值
[x,fval] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub);
figure();
subplot(2,2,[1,2])
yyaxis left
plot(1:14,x,'linewidth',1.5,'color','b');hold on;          % 最合理的压力规划
grid on;
ylabel('stress level');
xlabel('days');
yyaxis right
whr_c = zeros(1,14);
for i = 1:14
    for k = 1:i
        whr_c(i) = whr_c(i)+polyval(p,x(k));
    end
end
plot(1:14,whr_c*upr,'linewidth',1.5);hold on;              % 累计白发率
ylabel('cumulative white hair rate');
title('Optimization result');

% hair simulation VS
hair = zeros(fix(nh_l),fix(nh_l))+Initial;                % fix取整,设置头发初始灰度
hair = uint8(hair);
wi = 1;                                                   % white i,变白了的头发的索引
bi = 1;                                                   % black i,变黑了的头发的索引
wx = randi(fix(nh_l),1,fix(nh_l));
wy = randi(fix(nh_l),1,fix(nh_l));
for t = 1:1:14                                            % 对每一天
    xx = fix((nh_l-1)*rand(upr*nh,1)+1);                   % 求该天被更新的头发坐标集合，rand生成数组
    yy = fix((nh_l-1)*rand(upr*nh,1)+1);
    for i = 1:fix(upr*nh)
        i_whr = polyval(p,x(t));
        if i_whr > 0
            hair(xx(i),yy(i)) = hair(xx(i),yy(i))+255*i_whr;
            wx(wi) = xx(i);                                % 记录下变白了的头发
            wy(wi) = yy(i);
            wi = wi+1;
        else
            hair(wx(bi),wy(bi)) = ...
                hair(wx(bi),wy(bi))+255*i_whr;            % 此前变白的头发将优先恢复
            bi = bi+1;
        end
    end
end
subplot(2,2,3)
imshow(hair);
xlabel('optimal result');

hair = zeros(fix(nh_l),fix(nh_l))+Initial;                % fix取整,设置头发初始灰度
hair = uint8(hair);
wi = 1;                                                   % white i,变白了的头发的索引
bi = 1;                                                   % black i,变黑了的头发的索引
wx = randi(fix(nh_l),1,fix(nh_l));
wy = randi(fix(nh_l),1,fix(nh_l));
for t = 1:1:14                                            % 对每一天
    xx = fix((nh_l-1)*rand(upr*nh,1)+1);                   % 求该天被更新的头发坐标集合，rand生成数组
    yy = fix((nh_l-1)*rand(upr*nh,1)+1);
    for i = 1:fix(upr*nh)
        i_whr = (sl1+sl2)/2;
        if i_whr > 0
            hair(xx(i),yy(i)) = hair(xx(i),yy(i))+255*i_whr;
            wx(wi) = xx(i);                                % 记录下变白了的头发
            wy(wi) = yy(i);
            wi = wi+1;
        else
            hair(wx(bi),wy(bi)) = ...
                hair(wx(bi),wy(bi))+255*i_whr;            % 此前变白的头发将优先恢复
            bi = bi+1;
        end
    end
end
subplot(2,2,4)
imshow(hair);
xlabel('Non-optimal result result');