close all; clear;
%% 设置
fontsize = 20;
%%
datafile = 'data/US_electricity.xlsx';
sheetName = '2025.1.12';
T = readtable(datafile, 'Sheet', sheetName, 'VariableNamingRule', 'preserve');
satellite = char(T.platform);  % 提取卫星列
plume_quality = T.plume_quality;
wind_Speed = T.wind_speed;
emission_epa = T.emission_epa;
ratio_rgb_outlier = T.ratio_rgb_outlier;

%% 筛选
idx = true(length(satellite),1);
for i = 1:length(satellite)
    %     idx(i) = strcmp(satellite(i,1:2), 'ZY');  % 筛选高分/资源卫星
end

%% 颜色
color = [241, 64, 64; % 红
    204, 153, 0;  % 橙
    140, 180, 2; % 浅绿
    51, 144, 2]/255; % 深绿
sizes = [1 1 1 1 1] * 25;
colorData = zeros(size(plume_quality, 1), 3);
for i = 1:size(plume_quality, 1)
    if(plume_quality(i)==0)
        colorData(i, :) = color(1, :);
    elseif(plume_quality(i)==0.5)
        colorData(i, :) = color(2, :);
    elseif(plume_quality(i)==0.7)
        colorData(i, :) = color(3, :);
    elseif(plume_quality(i)==1)
        colorData(i, :) = color(4, :);
    else
        fprintf("Error input！");
    end
end
%% 风速与排放量
X1 = [wind_Speed(idx), emission_epa(idx)];
Y1 = classifyEmission(plume_quality(idx));
weights = adjustWeights(plume_quality(idx), 1);
% 调节模型
k1 = 30;
b1 = 220;
temp = X1(:, 2) - (k1 * X1(:, 1) + b1);
predLabels1 = false(size(temp,1), 1);
predLabels1(temp >= 0) = false;
predLabels1(temp < 0) = true;
score = ClassifyAssessment(Y1, predLabels1);

% 绘制风速与排放量的分类图
f1 = figure('position', [20 100 800 700]);
hold on;

XLIM1 = [0 11]; % 风速范围
YLIM2 = [0 2600]; % 排放量范围

% 绘制决策边界
xrange1 = linspace(XLIM1(1), XLIM1(2), 100);
% k1 = -line1.coef(1) / line1.coef(2);
% b1 = -line1.intercept / line1.coef(2);
yrange1 = k1 * xrange1 + b1;  % 使用方程 y = -(w1*x + b) / w2

% 在决策边界上方绘制浅绿色区域
h1 = fill([xrange1, fliplr(xrange1)], [yrange1, YLIM2(2) * ones(1, length(xrange1))], ...
    [204/255, 234/255, 141/255], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
h1.DisplayName = '';

% 在决策边界下方绘制浅色区域
h2 = fill([fliplr(xrange1), xrange1], [YLIM2(1) * ones(1, length(xrange1)), yrange1], ...
    [242/255, 227/255, 213/255], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
h2.DisplayName = '';

plot(xrange1, yrange1, 'k--', 'LineWidth', 1.5); % 绘制决策边界

h = gscatter(X1(:, 1), X1(:, 2), plume_quality(idx), color, '.', sizes); % 绘制散点

xlabel('Wind Speed (m/s)');
ylabel('EPA Emission Rate (ton/h)');

% 图例
% lgd = legend({'Unquantifiable', 'No Emission', 'Low Confidence',...
%     'Medium Confidence', 'High Confidence', 'Detect Limitation'});
lgd = legend(h, {'NE', 'LC', 'MC', 'HC'});

lgd.Position = [0.7, 0.7, 0.2, 0.2];  % 根据需要调整此位置，左下角为 (0,0)
lgd.FontSize = 12; % 设置字体
lgd.Box = 'off'; % 去掉图例背景框

% 在图窗中添加文本框，显示模型相关参数
text(0.65, 0.45, sprintf('y = %.1fx + %.0f', k1, b1), 'Units', 'normalized', 'FontSize', fontsize);
text(0.65, 0.35, sprintf('Accuracy: %.2f', score.accuracy), 'Units', 'normalized', 'FontSize', fontsize);
text(0.65, 0.25, sprintf('F1 score: %.2f', score.f1Score), 'Units', 'normalized', 'FontSize', fontsize);
% text(0.85, 0.55, sprintf('Pre: %.2f', score.precision), 'Units', 'normalized', 'FontSize', 12);
% text(0.85, 0.5, sprintf('Rec: %.2f', score.recall), 'Units', 'normalized', 'FontSize', 12);
box on;
axis square;
hold off;
xlim(XLIM1); % 风速范围
ylim(YLIM2); % 排放量范围
% 设置所有字体为 Calibri
set(gca, 'FontName', 'Calibri', 'FontSize', fontsize);
set(gcf, 'DefaultTextFontName', 'Calibri', 'DefaultTextFontSize', fontsize);
set(gcf, 'DefaultLegendFontName', 'Calibri', 'DefaultLegendFontSize', fontsize);

%% 背景与排放量
X2 = [ratio_rgb_outlier(idx), emission_epa(idx)];
Y2 = classifyEmission(plume_quality(idx));
weights = adjustWeights(plume_quality(idx), 2);

% 调节模型
k2 = 3200;
b2 = 230;
temp = X2(:, 2) - (k2 * X2(:, 1) + b2);
predLabels2 = false(size(temp,1), 1);
predLabels2(temp >= 0) = false;
predLabels2(temp < 0) = true;
score = ClassifyAssessment(Y2, predLabels2);

% 绘制风速与排放量的分类图
f2 = figure('position', [520 100 800 700]);
hold on;

XLIM2 = [0 0.1]; % 背景范围
YLIM2 = [0 2600]; % 排放量范围

% 绘制决策边界
xrange2 = linspace(XLIM2(1), XLIM2(2), 100);
yrange2 = k2 * xrange2 + b2;  % 使用方程 y = -(w1*x + b) / w2

% 在决策边界上方绘制浅绿色区域
h1 = fill([xrange2, fliplr(xrange2)], [yrange2, YLIM2(2) * ones(1, length(xrange2))], ...
    [204/255, 234/255, 141/255], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
h1.DisplayName = '';

% 在决策边界下方绘制浅色区域
h2 = fill([fliplr(xrange2), xrange2], [YLIM2(1) * ones(1, length(xrange2)), yrange2], ...
    [242/255, 227/255, 213/255], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
h2.DisplayName = '';

plot(xrange2, yrange2, 'k--', 'LineWidth', 1.5); % 绘制决策边界

h = gscatter(X2(:, 1), X2(:, 2), plume_quality(idx), color, '.', sizes); % 绘制散点

xlabel('Background Outlier Rate');
ylabel('EPA Emission Rate (ton/h)');

% 图例
lgd = legend(h, {'NE', 'LC','MC', 'HC'});

lgd.Position = [0.7, 0.7, 0.2, 0.2];  % 根据需要调整此位置，左下角为 (0,0)
lgd.FontSize = 12; % 设置字体
lgd.Box = 'off'; % 去掉图例背景框

% 在图窗中添加文本框，显示模型相关参数
text(0.65, 0.45, sprintf('y = %.0fx + %.0f', k2, b2), 'Units', 'normalized', 'FontSize', fontsize);
text(0.65, 0.35, sprintf('Accuracy: %.2f', score.accuracy), 'Units', 'normalized', 'FontSize', fontsize);
text(0.65, 0.25, sprintf('F1 score: %.2f', score.f1Score), 'Units', 'normalized', 'FontSize', fontsize);
% text(0.85, 0.55, sprintf('Pre: %.2f', score.precision), 'Units', 'normalized', 'FontSize', 12);
% text(0.85, 0.5, sprintf('Rec: %.2f', score.recall), 'Units', 'normalized', 'FontSize', 12);
box on;
axis square;
hold off;
xlim(XLIM2); % 背景范围
ylim(YLIM2); % 排放量范围
% 设置所有字体为 Calibri，字体大小为 12
set(gca, 'FontName', 'Calibri', 'FontSize', fontsize);
set(gcf, 'DefaultTextFontName', 'Calibri', 'DefaultTextFontSize', fontsize);
set(gcf, 'DefaultLegendFontName', 'Calibri', 'DefaultLegendFontSize', fontsize);

%% 三维散点图
f3 = figure('position', [820 100 800 700]);
hold on;

% 绘制三维散点图
scatter3(wind_Speed(idx), ratio_rgb_outlier(idx), emission_epa(idx), 50, colorData(idx, :), 'filled');
X3 = [wind_Speed(idx), ratio_rgb_outlier(idx), emission_epa(idx)];
Y3 = classifyEmission(plume_quality(idx));

% 定义 x1 和 x2 的范围
x1 = linspace(XLIM1(1), XLIM1(2), 50);  % 可根据需要调整范围
x2 = linspace(XLIM2(1), XLIM2(2), 50);  % 可根据需要调整范围

% 生成 x1 和 x2 的网格
[X1, X2] = meshgrid(x1, x2);

ka = 30;
kb = 3200;
b = 134;

% 计算平面 y 值
Y = ka * X1 + kb * X2 + b;
temp = X3(:, 3) - (ka * X3(:, 1) + kb * X3(:, 2) + b);
predLabels3 = false(size(temp,1), 1);
predLabels3(temp >= 0) = false;
predLabels3(temp < 0) = true;
score = ClassifyAssessment(Y3, predLabels3);

% 绘制半透明平面并添加光照效果
surfHandle = surf(X1, X2, Y);
surfHandle.FaceAlpha = 0.8;        % 设置透明度 (0.0 为全透明, 1.0 为不透明)
surfHandle.EdgeColor = 'none';     % 隐藏网格线
lighting gouraud;                  % 设置光照类型
camlight('headlight');             % 添加灯光，使平面更有立体感
% 设置坐标轴标题
xlabel('Wind Speed (m/s)');
ylabel('Outlier Rate');
zlabel('EPA Emission Rate (t/h)');

colormap('gray')
colorbar('off');
% 在图窗中添加文本框，显示模型相关参数
text(0.48, 0.40, sprintf('y = %.0f*x_1 + %.0f*x_2 + %.0f', ka, kb, b), 'Units', 'normalized', 'FontSize', fontsize);
text(0.70, 0.57, sprintf('Accuracy: %.2f', score.accuracy), 'Units', 'normalized', 'FontSize', fontsize);
text(0.70, 0.5, sprintf('F1 score: %.2f', score.f1Score), 'Units', 'normalized', 'FontSize', fontsize);
xticks(0:2:11); % 风速范围

% 设置所有字体为 Calibri
set(gca, 'FontName', 'Calibri', 'FontSize', fontsize);
set(gcf, 'DefaultTextFontName', 'Calibri', 'DefaultTextFontSize', fontsize);
set(gcf, 'DefaultLegendFontName', 'Calibri', 'DefaultLegendFontSize', fontsize);

% 设置视角
view(30, 7);
grid on;
box on;
hold off;
%% 绘制饼图
f4 = figure('position', [1120 100 800 700]);

% 定义质量等级和对应的标签
quality_levels = [1, 0.7, 0.5, 0];
labels = {'High Confidence', 'Medium Confidence', 'Low Confidence', 'No Emission'};

% 初始化计数器
counts = zeros(size(quality_levels));

% 统计每个质量等级的数量
for i = 1:length(quality_levels)
    counts(i) = sum(plume_quality == quality_levels(i));
end
percentage = counts/sum(counts)*100;
percentage_str = cell(length(quality_levels), 1);
for i = 1:length(quality_levels)
percentage_str{i} = char(sprintf("%.1f%%", percentage(i)));
end
% 定义红色系颜色
colors = {
    '#339002';
    '#8CB402';
    '#CC9900';
    '#F14040'};

% 绘制饼图并获取句柄
h = pie(counts, percentage_str);

% 设置颜色
for i = 1:length(h)/2
    set(h(2*i-1), 'FaceColor', colors{i}); % 设置饼片颜色
    set(h(2*i-1), 'EdgeColor', 'k');       % 设置饼片边缘颜色为黑色
    set(h(2*i-1), 'LineWidth', 1.5);         % 设置饼片边缘宽度
end

% 添加数据标签
for i = 1:2:length(h)
    % 在饼图上添加文本标签
    text_handle = findobj(f4, 'Type', 'text', 'Tag', '');
    if ~isempty(text_handle)
        set(text_handle((i+1)/2), 'Color', 'k', 'FontSize', fontsize);
    end
end

% 显示图例
lgd = legend(labels, 'Location', 'southoutside'); % 把图例放置在南部
lgd.NumColumns = 2; % 设置图例有两列
lgd.Box = 'off'; % 移除图例外边框
lgd.Position = [0.25 0.05 0.6112 0.0864];
set(gca, 'FontName', 'Arial', 'FontSize', fontsize);
%% 导出图像
print(f1,'C:\Users\huayi\Desktop\NC论文\Figure5\g1','-dpng','-r600');
print(f2,'C:\Users\huayi\Desktop\NC论文\Figure5\g2','-dpng','-r600');
print(f3,'C:\Users\huayi\Desktop\NC论文\Figure5\g3','-dpng','-r600');
print(f4,'C:\Users\huayi\Desktop\NC论文\Figure5\g4','-dpng','-r600');

%%
function [predLabels, line] = LinearClassifier(X, Y, weights) % 线性分类模型
% 训练线性分类器（线性判别分析）
ldaModel = fitcdiscr(X, Y, 'DiscrimType', 'linear', 'Weights', weights);
predLabels = predict(ldaModel, X);

% 获取线性分类器的系数
line.coef = ldaModel.Coeffs(1,2).Linear;  % 获取线性边界的系数
line.intercept = ldaModel.Coeffs(1,2).Const;  % 获取常数项
end

%%
function score = ClassifyAssessment(Y, predLabels) % 分类效果
% 计算混淆矩阵
cm = confusionmat(Y, predLabels);

% 提取 TP, TN, FP, FN
TP = cm(2, 2);  % 真正例
TN = cm(1, 1);  % 真反例
FP = cm(1, 2);  % 假正例
FN = cm(2, 1);  % 假反例

% 计算准确率
score.accuracy = (TP + TN) / (TP + TN + FP + FN);
disp(['Accuracy: ', num2str(score.accuracy)]);

% 计算精确率
score.precision = TP / (TP + FP);
disp(['Precision: ', num2str(score.precision)]);

% 计算召回率
score.recall = TP / (TP + FN);
disp(['Recall: ', num2str(score.recall)]);

% 计算 F1 分数
score.f1Score = 2 * (score.precision * score.recall) / (score.precision + score.recall);
disp(['F1 Score: ', num2str(score.f1Score)]);
end
%%
function Y = classifyEmission(Emission)
    % 根据 Emission 值分类
    Y = true(size(Emission));
    Y(Emission == 0.7) = false;
    Y(Emission == 1) = false;
    Y(Emission == -1) = false;
end
%%
function weights = adjustWeights(Emission, choice)
if(choice == 1)
    % 根据 Emission 值调整权重 Wind & EPA
    weights = ones(size(Emission));
    weights(Emission == 0) = 3;
    weights(Emission == 0.5) = 2;
    weights(Emission == 0.7) = 1;
    weights(Emission == 1) = 1;
elseif(choice == 2)
    % 根据 Emission 值调整权重 BKG & EPA 
    weights = ones(size(Emission));
    weights(Emission == 0) = 1;
    weights(Emission == 0.5) = 1;
    weights(Emission == 0.7) = 1;
    weights(Emission == 1) = 1;
end
end
