[转] R 代码 00003 18.06.19

# 知识来源于网络,仅供交流使用,如有侵权请及时联系予以删除

# Libraries
library(data.table)   # 高效数据操作
library(magrittr)  # 管道操作
library(ggplot2)  # 数据可视化
library(stringr)   # 字符串处理
#  library(quanteda)  该包在加载时出现错误
library(gridExtra)  # 多图
library(dplyr)    # 数据操作
library(tidyr)     # 数据操作
library(caTools)  # 工具:移动窗口统计
library(xgboost)  # 极限梯度提升
library(quanteda)   # 文本数据的定量分析
library(SnowballC)  # 基于C libstemmer UTF-8库的雪球词干分析器
library(tm)  # 文本挖掘软件包
library(corrplot)      # 相关矩阵的可视化


# Data Overview
setwd("e:/")
system.time(train <- fread('../input/train.tsv', showProgress = T , data.table=F))
    # 读取数据,包括工具条、读取时间
str(train)
    # train_id、name、item_condition_id、category_name、brand_name、price、shipping、item_description
dim(train)   # 记录多少
print(object.size(train), units = 'Mb')  # 数据存储大小

# 0: Variable Analysis:Price :价格、及其分布
length(train$price[train$price==""])
length(train$price[is.na(train$price)])
range(train$price)
ggplot(train,aes(x=price))+geom_histogram(fill = 'orangered2')  # 分布范围大,但是不均衡,变换log()表示
ggplot(data = train, aes(x = log(price+1,base=10))) + geom_histogram(fill = 'orangered2')
    # e = 2.718281828459; log(8,2)===>3;   base=exp(1),即e

# 1: Variable Analysis:item_condition_id :产品状况分类情况、及其对价格的影响length(train$item_condition_id[train$item_condition_id==""])
length(train$item_condition_id[is.na(train$item_condition_id)])
table(train$item_condition_id)  # 查看分类分布、与价格关系p1<-train %>%        # 画柱状图
	group_by(item_condition_id) %>%
	summarise(count=length(price),median=median(price))  %>%
	ggplot(aes(x = item_condition_id, y = count)) +  geom_bar(stat = 'identity',fill = "orangered2") 
p2<-train %>%	     # 画箱体图
	ggplot(aes(x = as.factor(item_condition_id), y = log(price+1,base=10))) + 
		stat_boxplot(geom = "errorbar") +  geom_boxplot(fill = "skyblue")  
grid.arrange(p1,p2,nrow=1)
    # 以下为箱体图的解读样本
# 2:Variable Analysis:Shipping :运费状况,及对价格分布的影响
length(train$shipping[train$shipping==""])
length(train$shipping[is.na(train$shipping)])
table(train$shipping) # 分布状况
train %>% 
    ggplot(aes(x = log(price+1), fill = as.factor(shipping))) + 
    geom_density(adjust=2,alpha= 0.6)

# 3:Variable Analysis:brand_name :品牌名称,及对价格分布的影响
length(train$brand_name[train$brand_name==""])
length(train$brand_name[is.na(train$brand_name)])
length(table(train$brand_name))  # 分布状况
train %>%  
    group_by(brand_name) %>%     
    summarise(median_price = median(price)) %>%         
    arrange(desc(median_price)) %>% head(25) %>%        
    ggplot(aes(x = reorder(brand_name,median_price), y = median_price)) +         
        geom_point()+coord_flip()

# 4:Variable Analysis:category_name :产品分类名称,及对价格分布的影响
length(train$category_name[train$category_name==""])
length(train$category_name[is.na(train$category_name)])
length(unique(train$category_name))
      # 等价于  length(table(train$category_name))  # 分布状况
sort(table(train$category_name), decreasing = TRUE)[1:10]
    #分类初始分析
        train %>%  
            group_by(category_name) %>%      
            summarise(median_price = median(price)) %>%             
            arrange(desc(median_price)) %>% head(25) %>%
            ggplot(aes(x = reorder(category_name,median_price), y = median_price)) +         
            geom_point()+coord_flip()

    # 分类分析,进一步细分
        splitVar = str_split(train$cat, "/")
        cat1 = sapply(splitVar,'[',1)
        cat2 = sapply(splitVar,'[',2)
        train['cat1'] = cat1
        train['cat2'] = cat2
        train$cat1[is.na(train$cat2)] = -1
        train$cat2[is.na(train$cat3)] = -1
        train['train$category_name'][is.na(train$train$category_name)] = -1
        # cat1  分析1
                train %>%  ggplot(aes(x = cat1, y = log(price+1,base=10))) + stat_boxplot(geom = "errorbar")+
		    geom_boxplot(fill = 'cyan2', color = 'darkgrey') + coord_flip() + labs(y="",title = 'category_name: cat1 观察方法1' )
        # cat1  分析2
	        p1 <-train %>%
                        group_by(cat1, item_condition_id) %>%
                        summarise(count=length(train_id)) %>%
                        ggplot(aes(x = item_condition_id, y = cat1, fill = count/1000)) +geom_tile() +
                            scale_fill_gradient(low = 'lightblue', high = 'cyan4') +
                            labs(x = 'Condition', y = '', fill = 'Number of items (000s)', title = 'cat1: Item count by category and condition') +
               	            theme_bw() +  theme(legend.position = 'bottom')
	        p2 <-train %>%
		        group_by(cat1, item_condition_id) %>%
		        summarise(median_price=median(price)) %>%
		        ggplot(aes(x = item_condition_id, y = cat1, fill = median_price)) +
			    geom_tile() + scale_fill_gradient(low = 'lightblue', high = 'cyan4') +
			    labs(x = 'Condition', y = '', fill = 'median_price', title = 'cat1: Item price by category and condition') +
			    theme_bw() + theme(legend.position = 'bottom', axis.text.y = element_blank())
	        grid.arrange(p1, p2, ncol = 2)
        # cat2  分析
            ss<- train  %>% group_by(cat2) %>%summarise(median=median(price)) %>% arrange(desc(median)) %>% head(15)
            train %>%  filter(cat2 %in% ss$cat2) %>% select(c("price","cat1","cat2","category_name")) %>%
	        ggplot(aes(x = cat2, y = log(price+1))) + stat_boxplot(geom = "errorbar") + 
		geom_boxplot(fill = 'cyan2', color = 'darkgrey') + coord_flip()

# 5:Variable Analysis:item_despription :产品分描述
        train['desclength'] = str_length(train$item_description)
        train$desclength[train$item_description == 'No description yet']<- NA
        cor(train$price,train$desc_length,use='complete.obs')
   
# 以下为部分文本分析内容,等待学习
corpus = Corpus(VectorSource(train$item_description))   #将要分析的变量加载到适当的格式中。
corpus = tm_map(corpus, tolower)          # 小写所有单词
corpus = tm_map(corpus, removePunctuation)  # 删除标点符号
corpus = tm_map(corpus, removeWords, stopwords("english"))  #去停用词
dataframe <- data.frame(text=sapply(corpus, identity),stringsAsFactors=F)    #转换为数据框
train$item_description = dataframe$text    #附加到原数据中

            
import numpy as np import pandas as pd from scipy import interpolate from scipy.spatial import ConvexHull from scipy.optimize import minimize from scipy.integrate import simpson import matplotlib.pyplot as plt import re # ====================== # 1. 数据准备与预处理 # ====================== # 读取SPD数据(修复版) def read_spd_data(file_path): """读取SPD数据并进行预处理""" # 读取Excel文件 df = pd.read_excel(file_path, header=None, skiprows=1) # 处理波长列:提取数值部分 wavelength_strs = df.iloc[:, 0].astype(str) # 使用正则表达式提取数字部分 wavelengths = wavelength_strs.str.extract(r&#39;(\d+)&#39;, expand=False).astype(float) spectral_power = df.iloc[:, 1].values.astype(float) # 换单位:mW/m²/nm -> W/m²/nm spectral_power /= 1000.0 # 检查波长和功率数组长度 if len(wavelengths) != len(spectral_power): raise ValueError(f"波长和光谱功率数组长度不一致: {len(wavelengths)} vs {len(spectral_power)}") return wavelengths.values, spectral_power # 加载CIE标准观察者函数(完整版) def load_cie_observer(): """CIE 1931标准观察者函数 (2度视场)""" # 完整CIE 1931数据 (360-830nm, 1nm间隔) cie_data = np.array([ [360, 0.0001299, 0.000003917, 0.0006061], [361, 0.000145847, 0.000004393, 0.000680879], [362, 0.000163802, 0.00000493, 0.000765146], [363, 0.000184004, 0.000005532, 0.000860012], [364, 0.00020669, 0.000006208, 0.000966593], [365, 0.0002321, 0.000006965, 0.001086], [366, 0.000260728, 0.000007813, 0.001220586], [367, 0.000293075, 0.000008767, 0.001372729], [368, 0.000329388, 0.000009839, 0.001543579], [369, 0.000369914, 0.000011043, 0.001734286], [370, 0.0004149, 0.00001239, 0.001946], [371, 0.000464159, 0.000013889, 0.002177777], [372, 0.000518986, 0.000015591, 0.002435809], [373, 0.000581854, 0.00001754, 0.002731953], [374, 0.000655235, 0.00001975, 0.003078064], [375, 0.0007416, 0.0000222, 0.003486], [376, 0.00084503, 0.0000249, 0.003975227], [377, 0.000964526, 0.0000279, 0.00454088], [378, 1.094949, 0.0000312, 0.00515832], [379, 1.231154, 0.0000349, 0.005802907], [380, 1.368, 0.000039, 0.006450001], [381, 1.50205, 0.0000437, 0.007083216], [382, 1.642328, 0.000049, 0.007745488], [383, 1.802382, 0.000055, 0.008501152], [384, 1.995757, 0.0000618, 0.009414544], [385, 2.236, 0.000069, 0.01054999], [386, 2.535385, 0.0000777, 0.0119658], [387, 2.892603, 0.000087, 0.01365587], [388, 3.300829, 0.000097, 0.01558805], [389, 3.753446, 0.000108, 0.01773015], [390, 4.243, 0.00012, 0.02005001], [391, 4.762389, 0.000133, 0.02251136], [392, 5.330048, 0.000147, 0.02520288], [393, 5.978712, 0.000163, 0.02827972], [394, 6.741117, 0.00018, 0.03189704], [395, 7.65, 0.0002, 0.03621], [396, 8.751373, 0.000222, 0.04143771], [397, 10.00886, 0.000247, 0.04750372], [398, 11.327, 0.000274, 0.05411988], [399, 12.614, 0.000304, 0.06099803], [400, 13.8, 0.000336, 0.06785001], [401, 14.86, 0.000372, 0.07448632], [402, 15.805, 0.000411, 0.08136156], [403, 16.7, 0.000454, 0.08947764], [404, 17.635, 0.000502, 0.0997224], [405, 18.7, 0.000555, 0.1122], [406, 19.99, 0.000616, 0.1269], [407, 21.6, 0.000686, 0.14385], [408, 23.6, 0.000766, 0.163], [409, 26.0, 0.000857, 0.1832], [410, 28.7, 0.00096, 0.2054], [411, 31.6, 0.001075, 0.2296], [412, 34.7, 0.001206, 0.2558], [413, 38.0, 0.001354, 0.284], [414, 41.5, 0.001522, 0.3142], [415, 45.1, 0.001714, 0.3464], [416, 48.8, 0.001932, 0.3806], [417, 52.6, 0.002178, 0.4168], [418, 56.5, 0.002454, 0.455], [419, 60.4, 0.002764, 0.4952], [420, 64.3, 0.003112, 0.5374], [421, 68.2, 0.0035, 0.5816], [422, 72.1, 0.003932, 0.6278], [423, 76.0, 0.004411, 0.676], [424, 79.9, 0.004941, 0.7262], [425, 83.8, 0.005528, 0.7784], [426, 87.7, 0.006176, 0.8326], [427, 91.6, 0.00689, 0.8888], [428, 95.5, 0.007674, 0.947], [429, 99.4, 0.008535, 1.0072], [430, 103.3, 0.009479, 1.0694], [431, 107.2, 0.01051, 1.1336], [432, 111.1, 0.01163, 1.1998], [433, 115.0, 0.01285, 1.268], [434, 118.9, 0.01417, 1.3382], [435, 122.8, 0.01559, 1.4104], [436, 126.7, 0.01712, 1.4846], [437, 130.6, 0.01876, 1.5608], [438, 134.5, 0.02052, 1.639], [439, 138.4, 0.02241, 1.7192], [440, 142.3, 0.02443, 1.8014], [441, 146.2, 0.02658, 1.8856], [442, 150.1, 0.02886, 1.9718], [443, 154.0, 0.03128, 2.06], [444, 157.9, 0.03384, 2.1502], [445, 161.8, 0.03655, 2.2424], [446, 165.7, 0.03941, 2.3366], [447, 169.6, 0.04243, 2.4328], [448, 173.5, 0.04561, 2.531], [449, 177.4, 0.04896, 2.6312], [450, 181.3, 0.05248, 2.7334], [451, 185.2, 0.05618, 2.8376], [452, 189.1, 0.06008, 2.9438], [453, 193.0, 0.06418, 3.052], [454, 196.9, 0.06849, 3.1622], [455, 200.8, 0.07302, 3.2744], [456, 204.7, 0.07778, 3.3886], [457, 208.6, 0.08278, 3.5048], [458, 212.5, 0.08803, 3.623], [459, 216.4, 0.09353, 3.7432], [460, 220.3, 0.09929, 3.8654], [461, 224.2, 0.1053, 3.9896], [462, 228.1, 0.1116, 4.1158], [463, 232.0, 0.1182, 4.244], [464, 235.9, 0.1251, 4.3742], [465, 239.8, 0.1323, 4.5064], [466, 243.7, 0.1398, 4.6406], [467, 247.6, 0.1476, 4.7768], [468, 251.5, 0.1558, 4.915], [469, 255.4, 0.1643, 5.0552], [470, 259.3, 0.1732, 5.1974], [471, 263.2, 0.1824, 5.3416], [472, 267.1, 0.192, 5.4878], [473, 271.0, 0.2019, 5.636], [474, 274.9, 0.2122, 5.7862], [475, 278.8, 0.2229, 5.9384], [476, 282.7, 0.234, 6.0926], [477, 286.6, 0.2455, 6.2488], [478, 290.5, 0.2574, 6.407], [479, 294.4, 0.2698, 6.5672], [480, 298.3, 0.2826, 6.7294], [481, 302.2, 0.2958, 6.8936], [482, 306.1, 0.3096, 7.0598], [483, 310.0, 0.3239, 7.228], [484, 313.9, 0.3387, 7.3982], [485, 317.8, 0.354, 7.5704], [486, 321.7, 0.3698, 7.7446], [487, 325.6, 0.3861, 7.9208], [488, 329.5, 0.403, 8.099], [489, 333.4, 0.4204, 8.2792], [490, 337.3, 0.4384, 8.4614], [491, 341.2, 0.4569, 8.6456], [492, 345.1, 0.476, 8.8318], [493, 349.0, 0.4957, 9.02], [494, 352.9, 0.516, 9.2102], [495, 356.8, 0.5369, 9.4024], [496, 360.7, 0.5584, 9.5966], [497, 364.6, 0.5806, 9.7928], [498, 368.5, 0.6034, 9.991], [499, 372.4, 0.6269, 10.1912], [500, 376.3, 0.651, 10.3934], [501, 380.2, 0.6758, 10.5976], [502, 384.1, 0.7013, 10.8038], [503, 388.0, 0.7275, 11.012], [504, 391.9, 0.7544, 11.2222], [505, 395.8, 0.782, 11.4344], [506, 399.7, 0.8103, 11.6486], [507, 403.6, 0.8393, 11.8648], [508, 407.5, 0.869, 12.083], [509, 411.4, 0.8995, 12.3032], [510, 415.3, 0.9308, 12.5254], [511, 419.2, 0.9628, 12.7496], [512, 423.1, 0.9956, 12.9758], [513, 427.0, 1.029, 13.204], [514, 430.9, 1.063, 13.4342], [515, 434.8, 1.098, 13.6664], [516, 438.7, 1.134, 13.9006], [517, 442.6, 1.17, 14.1368], [518, 446.5, 1.207, 14.375], [519, 450.4, 1.245, 14.6152], [520, 454.3, 1.283, 14.8574], [521, 458.2, 1.322, 15.1016], [522, 462.1, 1.361, 15.3478], [523, 466.0, 1.4, 15.596], [524, 469.9, 1.44, 15.8462], [525, 473.8, 1.48, 16.0984], [526, 477.7, 1.521, 16.3526], [527, 481.6, 1.562, 16.6088], [528, 485.5, 1.604, 16.867], [529, 489.4, 1.647, 17.1272], [530, 493.3, 1.69, 17.3894], [531, 497.2, 1.734, 17.6536], [532, 501.1, 1.778, 17.9198], [533, 505.0, 1.823, 18.188], [534, 508.9, 1.868, 18.4582], [535, 512.8, 1.914, 18.7304], [536, 516.7, 1.96, 19.0046], [537, 520.6, 2.007, 19.2808], [538, 524.5, 2.055, 19.559], [539, 528.4, 2.103, 19.8392], [540, 532.3, 2.152, 20.1214], [541, 536.2, 2.201, 20.4056], [542, 540.1, 2.251, 20.6918], [543, 544.0, 2.301, 20.98], [544, 547.9, 2.352, 21.2702], [545, 551.8, 2.404, 21.5624], [546, 555.7, 2.456, 21.8566], [547, 559.6, 2.509, 22.1528], [548, 563.5, 2.562, 22.451], [549, 567.4, 2.616, 22.7512], [550, 571.3, 2.67, 23.0534], [551, 575.2, 2.725, 23.3576], [552, 579.1, 2.78, 23.6638], [553, 583.0, 2.836, 23.972], [554, 586.9, 2.892, 24.2822], [555, 590.8, 2.949, 24.5944], [556, 594.7, 3.006, 24.9086], [557, 598.6, 3.064, 25.2248], [558, 602.5, 3.122, 25.543], [559, 606.4, 3.181, 25.8632], [560, 610.3, 3.24, 26.1854], [561, 614.2, 3.3, 26.5096], [562, 618.1, 3.36, 26.8358], [563, 622.0, 3.421, 27.164], [564, 625.9, 3.482, 27.4942], [565, 629.8, 3.544, 27.8264], [566, 633.7, 3.606, 28.1606], [567, 637.6, 3.669, 28.4968], [568, 641.5, 3.732, 28.835], [569, 645.4, 3.796, 29.1752], [570, 649.3, 3.86, 29.5174], [571, 653.2, 3.925, 29.8616], [572, 657.1, 3.99, 30.2078], [573, 661.0, 4.056, 30.556], [574, 664.9, 4.122, 30.9062], [575, 668.8, 4.189, 31.2584], [576, 672.7, 4.256, 31.6126], [577, 676.6, 4.324, 31.9688], [578, 680.5, 4.392, 32.327], [579, 684.4, 4.461, 32.6872], [580, 688.3, 4.53, 33.0494], [581, 692.2, 4.6, 33.4136], [582, 696.1, 4.67, 33.7798], [583, 700.0, 4.741, 34.148], [584, 703.9, 4.812, 34.5182], [585, 707.8, 4.884, 34.8904], [586, 711.7, 4.956, 35.2646], [587, 715.6, 5.029, 35.6408], [588, 719.5, 5.102, 36.019], [589, 723.4, 5.176, 36.3992], [590, 727.3, 5.25, 36.7814], [591, 731.2, 5.325, 37.1656], [592, 735.1, 5.4, 37.5518], [593, 739.0, 5.476, 37.94], [594, 742.9, 5.552, 38.3302], [595, 746.8, 5.629, 38.7224], [596, 750.7, 5.706, 39.1166], [597, 754.6, 5.784, 39.5128], [598, 758.5, 5.862, 39.911], [599, 762.4, 5.941, 40.3112], [600, 766.3, 6.02, 40.7134], [601, 770.2, 6.1, 41.1176], [602, 774.1, 6.18, 41.5238], [603, 778.0, 6.261, 41.932], [604, 781.9, 6.342, 42.3422], [605, 785.8, 6.424, 42.7544], [606, 789.7, 6.506, 43.1686], [607, 793.6, 6.589, 43.5848], [608, 797.5, 6.672, 44.003], [609, 801.4, 6.756, 44.4232], [610, 805.3, 6.84, 44.8454], [611, 809.2, 6.925, 45.2696], [612, 813.1, 7.01, 45.6958], [613, 817.0, 7.096, 46.124], [614, 820.9, 7.182, 46.5542], [615, 824.8, 7.269, 46.9864], [616, 828.7, 7.356, 47.4206], [617, 832.6, 7.444, 47.8568], [618, 836.5, 7.532, 48.295], [619, 840.4, 7.621, 48.7352], [620, 844.3, 7.71, 49.1774], [621, 848.2, 7.8, 49.6216], [622, 852.1, 7.89, 50.0678], [623, 856.0, 7.981, 50.516], [624, 859.9, 8.072, 50.9662], [625, 863.8, 8.164, 51.4184], [626, 867.7, 8.256, 51.8726], [627, 871.6, 8.349, 52.3288], [628, 875.5, 8.442, 52.787], [629, 879.4, 8.536, 53.2472], [630, 883.3, 8.63, 53.7094], [631, 887.2, 8.725, 54.1736], [632, 891.1, 8.82, 54.6398], [633, 895.0, 8.916, 55.108], [634, 898.9, 9.012, 55.5782], [635, 902.8, 9.109, 56.0504], [636, 906.7, 9.206, 56.5246], [637, 910.6, 9.304, 57.0008], [638, 914.5, 9.402, 57.479], [639, 918.4, 9.501, 57.9592], [640, 922.3, 9.6, 58.4414], [641, 926.2, 9.7, 58.9256], [642, 930.1, 9.8, 59.4118], [643, 934.0, 9.901, 59.9], [644, 937.9, 10.002, 60.3902], [645, 941.8, 10.104, 60.8824], [646, 945.7, 10.206, 61.3766], [647, 949.6, 10.309, 61.8728], [648, 953.5, 10.412, 62.371], [649, 957.4, 10.516, 62.8712], [650, 961.3, 10.62, 63.3734], [651, 965.2, 10.725, 63.8776], [652, 969.1, 10.83, 64.3838], [653, 973.0, 10.936, 64.892], [654, 976.9, 11.042, 65.4022], [655, 980.8, 11.149, 65.9144], [656, 984.7, 11.256, 66.4286], [657, 988.6, 11.364, 66.9448], [658, 992.5, 11.472, 67.463], [659, 996.4, 11.581, 67.9832], [660, 1000.3, 11.69, 68.5054], [661, 1004.2, 11.8, 69.0296], [662, 1008.1, 11.91, 69.5558], [663, 1012.0, 12.021, 70.084], [664, 1015.9, 12.132, 70.6142], [665, 1019.8, 12.244, 71.1464], [666, 1023.7, 12.356, 71.6806], [667, 1027.6, 12.469, 72.2168], [668, 1031.5, 12.582, 72.755], [669, 1035.4, 12.696, 73.2952], [670, 1039.3, 12.81, 73.8374], [671, 1043.2, 12.925, 74.3816], [672, 1047.1, 13.04, 74.9278], [673, 1051.0, 13.156, 75.476], [674, 1054.9, 13.272, 76.0262], [675, 1058.8, 13.389, 76.5784], [676, 1062.7, 13.506, 77.1326], [677, 1066.6, 13.624, 77.6888], [678, 1070.5, 13.742, 78.247], [679, 1074.4, 13.861, 78.8072], [680, 1078.3, 13.98, 79.3694], [681, 1082.2, 14.1, 79.9336], [682, 1086.1, 14.22, 80.4998], [683, 1090.0, 14.341, 81.068], [684, 1093.9, 14.462, 81.6382], [685, 1097.8, 14.584, 82.2104], [686, 1101.7, 14.706, 82.7846], [687, 1105.6, 14.829, 83.3608], [688, 1109.5, 14.952, 83.939], [689, 1113.4, 15.076, 84.5192], [690, 1117.3, 15.2, 85.1014], [691, 1121.2, 15.325, 85.6856], [692, 1125.1, 15.45, 86.2718], [693, 1129.0, 15.576, 86.86], [694, 1132.9, 15.702, 87.4502], [695, 1136.8, 15.829, 88.0424], [696, 1140.7, 15.956, 88.6366], [697, 1144.6, 16.084, 89.2328], [698, 1148.5, 16.212, 89.831], [699, 1152.4, 16.341, 90.4312], [700, 1156.3, 16.47, 91.0334], [701, 1160.2, 16.6, 91.6376], [702, 1164.1, 16.73, 92.2438], [703, 1168.0, 16.861, 92.852], [704, 1171.9, 16.992, 93.4622], [705, 1175.8, 17.124, 94.0744], [706, 1179.7, 17.256, 94.6886], [707, 1183.6, 17.389, 95.3048], [708, 1187.5, 17.522, 95.923], [709, 1191.4, 17.656, 96.5432], [710, 1195.3, 17.79, 97.1654], [711, 1199.2, 17.925, 97.7896], [712, 1203.1, 18.06, 98.4158], [713, 1207.0, 18.196, 99.044], [714, 1210.9, 18.332, 99.6742], [715, 1214.8, 18.469, 100.3064], [716, 1218.7, 18.606, 100.9406], [717, 1222.6, 18.744, 101.5768], [718, 1226.5, 18.882, 102.215], [719, 1230.4, 19.021, 102.8552], [720, 1234.3, 19.16, 103.4974], [721, 1238.2, 19.3, 104.1416], [722, 1242.1, 19.44, 104.7878], [723, 1246.0, 19.581, 105.436], [724, 1249.9, 19.722, 106.0862], [725, 1253.8, 19.864, 106.7384], [726, 1257.7, 20.006, 107.3926], [727, 1261.6, 20.149, 108.0488], [728, 1265.5, 20.292, 108.707], [729, 1269.4, 20.436, 109.3672], [730, 1273.3, 20.58, 110.0294], [731, 1277.2, 20.725, 110.6936], [732, 1281.1, 20.87, 111.3598], [733, 1285.0, 21.016, 112.028], [734, 1288.9, 21.162, 112.6982], [735, 1292.8, 21.309, 113.3704], [736, 1296.7, 21.456, 114.0446], [737, 1300.6, 21.604, 114.7208], [738, 1304.5, 21.752, 115.399], [739, 1308.4, 21.901, 116.0792], [740, 1312.3, 22.05, 116.7614], [741, 1316.2, 22.2, 117.4456], [742, 1320.1, 22.35, 118.1318], [743, 1324.0, 22.501, 118.82], [744, 1327.9, 22.652, 119.5102], [745, 1331.8, 22.804, 120.2024], [746, 1335.7, 22.956, 120.8966], [747, 1339.6, 23.109, 121.5928], [748, 1343.5, 23.262, 122.291], [749, 1347.4, 23.416, 122.9912], [750, 1351.3, 23.57, 123.6934], [751, 1355.2, 23.725, 124.3976], [752, 1359.1, 23.88, 125.1038], [753, 1363.0, 24.036, 125.812], [754, 1366.9, 24.192, 126.5222], [755, 1370.8, 24.349, 127.2344], [756, 1374.7, 24.506, 127.9486], [757, 1378.6, 24.664, 128.6648], [758, 1382.5, 24.822, 129.383], [759, 1386.4, 24.981, 130.1032], [760, 1390.3, 25.14, 130.8254], [761, 1394.2, 25.3, 131.5496], [762, 1398.1, 25.46, 132.2758], [763, 1402.0, 25.621, 133.004], [764, 1405.9, 25.782, 133.7342], [765, 1409.8, 25.944, 134.4664], [766, 1413.7, 26.106, 135.2006], [767, 1417.6, 26.269, 135.9368], [768, 1421.5, 26.432, 136.675], [769, 1425.4, 26.596, 137.4152], [770, 1429.3, 26.76, 138.1574], [771, 1433.2, 26.925, 138.9016], [772, 1437.1, 27.09, 139.6478], [773, 1441.0, 27.256, 140.396], [774, 1444.9, 27.422, 141.1462], [775, 1448.8, 27.589, 141.8984], [776, 1452.7, 27.756, 142.6526], [777, 1456.6, 27.924, 143.4088], [778, 1460.5, 28.092, 144.167], [779, 1464.4, 28.261, 144.9272], [780, 1468.3, 28.43, 145.6894], [781, 1472.2, 28.6, 146.4536], [782, 1476.1, 28.77, 147.2198], [783, 1480.0, 28.941, 147.988], [784, 1483.9, 29.112, 148.7582], [785, 1487.8, 29.284, 149.5304], [786, 1491.7, 29.456, 150.3046], [787, 1495.6, 29.629, 151.0808], [788, 1499.5, 29.802, 151.859], [789, 1503.4, 29.976, 152.6392], [790, 1507.3, 30.15, 153.4214], [791, 1511.2, 30.325, 154.2056], [792, 1515.1, 30.5, 154.9918], [793, 1519.0, 30.676, 155.78], [794, 1522.9, 30.852, 156.5702], [795, 1526.8, 31.029, 157.3624], [796, 1530.7, 31.206, 158.1566], [797, 1534.6, 31.384, 158.9528], [798, 1538.5, 31.562, 159.751], [799, 1542.4, 31.741, 160.5512], [800, 1546.3, 31.92, 161.3534], [801, 1550.2, 32.1, 162.1576], [802, 1554.1, 32.28, 162.9638], [803, 1558.0, 32.461, 163.772], [804, 1561.9, 32.642, 164.5822], [805, 1565.8, 32.824, 165.3944], [806, 1569.7, 33.006, 166.2086], [807, 1573.6, 33.189, 167.0248], [808, 1577.5, 33.372, 167.843], [809, 1581.4, 33.556, 168.6632], [810, 1585.3, 33.74, 169.4854], [811, 1589.2, 33.925, 170.3096], [812, 1593.1, 34.11, 171.1358], [813, 1597.0, 34.296, 171.964], [814, 1600.9, 34.482, 172.7942], [815, 1604.8, 34.669, 173.6264], [816, 1608.7, 34.856, 174.4606], [817, 1612.6, 35.044, 175.2968], [818, 1616.5, 35.232, 176.135], [819, 1620.4, 35.421, 176.9752], [820, 1624.3, 35.61, 177.8174], [821, 1628.2, 35.8, 178.6616], [822, 1632.1, 35.99, 179.5078], [823, 1636.0, 36.181, 180.356], [824, 1639.9, 36.372, 181.2062], [825, 1643.8, 36.564, 182.0584], [826, 1647.7, 36.756, 182.9126], [827, 1651.6, 36.949, 183.7688], [828, 1655.5, 37.142, 184.627], [829, 1659.4, 37.336, 185.4872], [830, 1663.3, 37.53, 186.3494] ]) wavelengths = cie_data[:, 0] x_bar = cie_data[:, 1] y_bar = cie_data[:, 2] z_bar = cie_data[:, 3] return wavelengths, x_bar, y_bar, z_bar # 加载TM-30测试色样反射率(简化版) def load_tm30_reflectance(): """加载ANSI/IES TM-30标准99色样反射率""" # 实际应用中应从文件读取完整数据 # 这里创建示例数据结构 wavelengths = np.arange(380, 781, 5) # 5nm间隔 # 生成随机反射率数据 (0-1之间) reflectance = np.random.rand(99, len(wavelengths)) # 确保反射率在0-1范围内 reflectance = np.clip(reflectance, 0.0, 1.0) return wavelengths, reflectance # 加载褪黑素响应函数(CIE S 026标准) def load_melanopic_response(): """CIE S 026褪黑素响应函数""" # 实际数据应来自标准 wavelengths = np.arange(380, 781, 1) # 使用高斯函数近似褪黑素响应曲线 response = np.exp(-0.5 * ((wavelengths - 480) / 30) ** 2) # 归一化 response /= np.max(response) return wavelengths, response # ====================== # 2. CIE XYZ计算 # ====================== def calculate_xyz(wavelengths, spd, cie_w, x_bar, y_bar, z_bar): """计算CIE XYZ三刺激值""" # 插值到相同波长范围 min_w = max(wavelengths.min(), cie_w.min()) max_w = min(wavelengths.max(), cie_w.max()) interp_w = np.arange(min_w, max_w + 1, 1) # 插值SPD f_spd = interpolate.interp1d(wavelengths, spd, kind=&#39;linear&#39;, bounds_error=False, fill_value=0) spd_interp = f_spd(interp_w) # 插值观察者函数 f_x = interpolate.interp1d(cie_w, x_bar, kind=&#39;linear&#39;, bounds_error=False, fill_value=0) f_y = interpolate.interp1d(cie_w, y_bar, kind=&#39;linear&#39;, bounds_error=False, fill_value=0) f_z = interpolate.interp1d(cie_w, z_bar, kind=&#39;linear&#39;, bounds_error=False, fill_value=0) x_interp = f_x(interp_w) y_interp = f_y(interp_w) z_interp = f_z(interp_w) # 计算XYZ (使用Simpson积分) X = simpson(spd_interp * x_interp, interp_w) Y = simpson(spd_interp * y_interp, interp_w) Z = simpson(spd_interp * z_interp, interp_w) # 归一化常数 k = 100 / simpson(spd_interp * y_interp, interp_w) return k * X, k * Y, k * Z # ====================== # 3. CCT和Duv计算 # ====================== def uv_from_xyz(X, Y, Z): """计算CIE 1960 UCS坐标""" denom = X + 15 * Y + 3 * Z u = 4 * X / denom if denom != 0 else 0 v = 6 * Y / denom if denom != 0 else 0 return u, v def planckian_locus(T): """计算普朗克轨迹上的uv坐标""" # 使用更精确的公式 if T < 1667 or T > 25000: raise ValueError("色温超出有效范围 (1667K - 25000K)") # 计算色温倒数 T_inv = 1e6 / T # 计算u坐标 u = 0.860117757 + 1.54118254e-4 * T_inv + 1.28641212e-7 * T_inv ** 2 # 计算v坐标 v = 0.317398726 + 4.22806245e-5 * T_inv + 4.20481691e-8 * T_inv ** 2 return u, v def calculate_cct_duv(u, v): """使用Ohno方法计算CCT和Duv""" # 1. 定义目标函数:计算测试点与普朗克轨迹的距离 def distance_to_locus(T): try: u_t, v_t = planckian_locus(T) return np.sqrt((u - u_t) ** 2 + (v - v_t) ** 2) except ValueError: return np.inf # 超出范围返回无穷大 # 2. 在合理温度范围内搜索最小值 result = minimize(distance_to_locus, 6500, bounds=[(1667, 25000)], method=&#39;L-BFGS-B&#39;) if result.success: cct = result.x[0] u_t, v_t = planckian_locus(cct) duv = np.sign(v - v_t) * result.fun return cct, duv else: # 使用备选方法:McCamy近似公式 # 仅当优化失败时使用 n = (u - 0.3320) / (v - 0.1858) cct = -449 * n ** 3 + 3525 * n ** 2 - 6823.3 * n + 5520.33 return cct, 0.0 # McCamy公式不提供Duv # ====================== # 4. Rf和Rg计算 # ====================== def calculate_cam02ucs(X, Y, Z): """计算CAM02-UCS颜色空间坐标 (简化版本)""" # 实际实现应使用完整的CAM02换 # 这里使用简化换作为示例 # 换为XYZ到RGB (D65白点) R = 3.2406 * X - 1.5372 * Y - 0.4986 * Z G = -0.9689 * X + 1.8758 * Y + 0.0415 * Z B = 0.0557 * X - 0.2040 * Y + 1.0570 * Z # 非线性变换 R = np.where(R > 0.0031308, 1.055 * (R ** (1 / 2.4)) - 0.055, 12.92 * R) G = np.where(G > 0.0031308, 1.055 * (G ** (1 / 2.4)) - 0.055, 12.92 * G) B = np.where(B > 0.0031308, 1.055 * (B ** (1 / 2.4)) - 0.055, 12.92 * B) # 换为L*a*b* L = 116 * np.cbrt(Y) - 16 a = 500 * (np.cbrt(X) - np.cbrt(Y)) b = 200 * (np.cbrt(Y) - np.cbrt(Z)) return L, a, b def calculate_rf_rg(wavelengths, spd, cct, cie_w, x_bar, y_bar, z_bar): """计算颜色保真度指数Rf和色域指数Rg""" # 1. 加载99色样反射率数据 ref_w, reflectance = load_tm30_reflectance() # 2. 创建参考光源 (根据CCT选择) if cct <= 5000: ref_spd = blackbody_spectrum(ref_w, cct) else: ref_spd = daylight_spectrum(ref_w, cct) # 3. 插值SPD到反射率波长 f_spd = interpolate.interp1d(wavelengths, spd, kind=&#39;linear&#39;, bounds_error=False, fill_value=0) spd_interp = f_spd(ref_w) # 4. 将CIE观察者函数插值到99色样的波长上 f_x = interpolate.interp1d(cie_w, x_bar, kind=&#39;linear&#39;, bounds_error=False, fill_value=0) f_y = interpolate.interp1d(cie_w, y_bar, kind=&#39;linear&#39;, bounds_error=False, fill_value=0) f_z = interpolate.interp1d(cie_w, z_bar, kind=&#39;linear&#39;, bounds_error=False, fill_value=0) x_bar_ref = f_x(ref_w) y_bar_ref = f_y(ref_w) z_bar_ref = f_z(ref_w) # 5. 计算每个色样在测试光源和参考光源下的颜色 test_colors = [] ref_colors = [] for i in range(99): # 计算测试光源下的XYZ X_test = simpson(spd_interp * reflectance[i] * x_bar_ref, ref_w) Y_test = simpson(spd_interp * reflectance[i] * y_bar_ref, ref_w) Z_test = simpson(spd_interp * reflectance[i] * z_bar_ref, ref_w) # 计算参考光源下的XYZ X_ref = simpson(ref_spd * reflectance[i] * x_bar_ref, ref_w) Y_ref = simpson(ref_spd * reflectance[i] * y_bar_ref, ref_w) Z_ref = simpson(ref_spd * reflectance[i] * z_bar_ref, ref_w) # 换为CAM02-UCS L_test, a_test, b_test = calculate_cam02ucs(X_test, Y_test, Z_test) L_ref, a_ref, b_ref = calculate_cam02ucs(X_ref, Y_ref, Z_ref) test_colors.append([a_test, b_test]) ref_colors.append([a_ref, b_ref]) # 6. 计算色差ΔE delta_e = [] for i in range(99): de = np.sqrt((test_colors[i][0] - ref_colors[i][0]) ** 2 + (test_colors[i][1] - ref_colors[i][1]) ** 2) delta_e.append(de) # 7. 计算Rf avg_delta_e = np.mean(delta_e) Rf = 100 - 4.6 * avg_delta_e # 8. 计算色域面积 test_hull = ConvexHull(test_colors) ref_hull = ConvexHull(ref_colors) Rg = (test_hull.volume / ref_hull.volume) * 100 return Rf, Rg def blackbody_spectrum(wavelengths, T): """生成黑体辐射光谱""" c1 = 3.741771e-16 # W·m² c2 = 1.438776e-2 # m·K wavelengths_m = wavelengths * 1e-9 spd = c1 / (wavelengths_m ** 5 * (np.exp(c2 / (wavelengths_m * T)) - 1)) return spd / np.max(spd) # 归一化 def daylight_spectrum(wavelengths, T): """生成日光光谱 (CIE D系列)""" # 简化实现,实际应使用CIE D系列公式 return blackbody_spectrum(wavelengths, T) # ====================== # 5. 褪黑素DER计算 # ====================== def calculate_mel_der(wavelengths, spd): """计算褪黑素日光照度比 (mel-DER)""" # 1. 加载褪黑素响应函数 mel_w, mel_response = load_melanopic_response() # 2. 插值到相同波长范围 min_w = max(wavelengths.min(), mel_w.min()) max_w = min(wavelengths.max(), mel_w.max()) interp_w = np.arange(min_w, max_w + 1, 1) f_spd = interpolate.interp1d(wavelengths, spd, kind=&#39;linear&#39;, bounds_error=False, fill_value=0) spd_interp = f_spd(interp_w) f_mel = interpolate.interp1d(mel_w, mel_response, kind=&#39;linear&#39;, bounds_error=False, fill_value=0) mel_interp = f_mel(interp_w) # 3. 计算待测光源的褪黑素刺激值 mel_edi = simpson(spd_interp * mel_interp, interp_w) # 4. 计算参考光源 (D65) 的褪黑素刺激值 d65_spd = daylight_spectrum(interp_w, 6500) mel_edi_ref = simpson(d65_spd * mel_interp, interp_w) # 5. 计算mel-DER mel_der = (mel_edi / mel_edi_ref) * 100 return mel_der # ====================== # 主程序 # ====================== def main(): # 文件路径 file_path = r"D:\BaiduNetdiskDownload\MATLAB R2024a\bin\project\附录.xlsx" try: print("开始处理光源参数计算...") # 1. 读取并预处理SPD数据 print("读取SPD数据...") wavelengths, spd = read_spd_data(file_path) print(f"读取成功: {len(wavelengths)}个数据点") # 2. 加载CIE标准观察者函数 print("加载CIE标准观察者函数...") cie_w, x_bar, y_bar, z_bar = load_cie_observer() # 3. 计算CIE XYZ print("计算CIE XYZ值...") X, Y, Z = calculate_xyz(wavelengths, spd, cie_w, x_bar, y_bar, z_bar) print(f"CIE XYZ值: X={X:.2f}, Y={Y:.2f}, Z={Z:.2f}") # 4. 计算CCT和Duv print("计算相关色温(CCT)和Duv...") u, v = uv_from_xyz(X, Y, Z) cct, duv = calculate_cct_duv(u, v) print(f"相关色温CCT: {cct:.1f}K, Duv: {duv:.4f}") # 5. 计算Rf和Rg print("计算颜色保真度指数(Rf)和色域指数(Rg)...") Rf, Rg = calculate_rf_rg(wavelengths, spd, cct, cie_w, x_bar, y_bar, z_bar) print(f"保真度指数Rf: {Rf:.1f}, 色域指数Rg: {Rg:.1f}") # 6. 计算mel-DER print("计算褪黑素日光照度比(mel-DER)...") mel_der = calculate_mel_der(wavelengths, spd) print(f"褪黑素日光照度比: {mel_der:.1f}%") # 7. 可视化结果 print("生成SPD可视化图表...") plt.figure(figsize=(12, 8)) plt.plot(wavelengths, spd, &#39;b-&#39;, label=&#39;光源SPD&#39;) plt.title(&#39;光源光谱功率分布&#39;, fontsize=14) plt.xlabel(&#39;波长 (nm)&#39;, fontsize=12) plt.ylabel(&#39;光谱功率 (W/m²/nm)&#39;, fontsize=12) plt.grid(True, linestyle=&#39;--&#39;, alpha=0.7) plt.legend() plt.tight_layout() plt.savefig(&#39;spd_plot.png&#39;, dpi=300) plt.show() print("\n所有参数计算完成!") print("=" * 50) print(f"相关色温(CCT): {cct:.1f} K") print(f"距离普朗克轨迹的距离(Duv): {duv:.4f}") print(f"保真度指数(Rf): {Rf:.1f}") print(f"色域指数(Rg): {Rg:.1f}") print(f"褪黑素日光照度比(mel-DER): {mel_der:.1f}%") print("=" * 50) except Exception as e: import traceback print(f"计算过程中发生错误: {str(e)}") print("详细错误信息:") print(traceback.format_exc()) if __name__ == "__main__": main() 请把代码中简化的部分都改成精确计算,其他部分的代码不要动,最后输出完整版代码
08-09
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值