diff --git a/afc/radiance/forecast2.py b/afc/radiance/forecast2.py deleted file mode 100644 index 144457d..0000000 --- a/afc/radiance/forecast2.py +++ /dev/null @@ -1,646 +0,0 @@ -# Advanced Fenestration Controller (AFC) Copyright (c) 2023, The -# Regents of the University of California, through Lawrence Berkeley -# National Laboratory (subject to receipt of any required approvals -# from the U.S. Dept. of Energy). All rights reserved. - -""""Advanced Fenestration Controller -Radiance forecasting module. -""" - -# pylint: disable=too-many-locals, too-many-instance-attributes, too-many-arguments -# pylint: disable=redefined-outer-name, invalid-name, too-many-statements -# pylint: disable=consider-using-dict-items, protected-access - -import os -import time -import multiprocessing as mp -from configparser import ConfigParser - -import numpy as np -import pandas as pd -import matplotlib.pyplot as plt -from frads import room, utils, matrix, methods - -from .maps import shade_map_0x6, shade_map_0x4 - -root = os.path.dirname(os.path.realpath(__file__)) - -mult = np.linalg.multi_dot - - -KFOMG = ( - [0.0238639257641843] - + [0.02332285973188108] * 8 - + [0.021916319185673504] * 16 - + [0.023622221568953158] * 20 - + [0.022385166034409335] * 24 - + [0.022385166034409376] * 24 - + [0.019685184640794287] * 24 - + [0.021916319185673487] * 16 - + [0.0175372336349361] * 12 -) - - -class Forecast: - """Radiance forecaster class for the AFC.""" - - def __init__( - self, - cfg_path, - regenerate=None, - facade_type="ec", - wpi_plot=False, - wpi_loc="23back", - location=None, - filestruct=None, - _test_difference=False, - dimensions=None, - ): - self.root = os.path.dirname(os.path.abspath(__file__)) - self.parse_config(cfg_path) - if regenerate: - self.remake_matrices = True - if filestruct: - self.filestruct = filestruct - self.bsdfd = self.filestruct["resources"] - self.matricesd = self.filestruct["matrices"] - if location: - self.lat = location["latitude"] - self.lon = location["longitude"] - self.timezone = location["timezone"] - self.elevation = location["elevation"] - self.orient = location["orient"] - if dimensions: - self.dims["width"] = dimensions["width"] - self.dims["depth"] = dimensions["depth"] - self.dims["height"] = dimensions["height"] - self.dims["window1"] = dimensions["window1"] - self.dims["window2"] = dimensions["window2"] - self.dims["window3"] = dimensions["window3"] - self.facade_type = facade_type - self.wpi_plot = wpi_plot - self.wpi_loc = wpi_loc - self._test_difference = _test_difference - self.get_paths() - self.make_room() - self.get_workflow() - self.make_matrices() - self.load_mtx() - self.wwr = float(cfg_path.split("room")[1].split("WWR")[0]) - if "shade" in self.facade_type and float(self.wwr) == 0.6: - self.shade_map = shade_map_0x6 - elif "shade" in self.facade_type and float(self.wwr) == 0.4: - self.shade_map = shade_map_0x4 - self.new_map = {} - self.init = True - - # Initialize variables - self.wpi_data = None - self.map_center = None - self.map_23 = None - self.n_w = None - self.n_d = None - self.date_time = None - - def parse_config(self, cfg_path): - """Function to parse configuration file.""" - - config = ConfigParser() - if not config.read(cfg_path): - raise ValueError( - 'The location of the "config" file is wrong. ' - + f"Check location at {cfg_path}." - ) - cfg = config._sections - site_info = cfg["Site"] - self.lat = site_info["latitude"] - self.lon = site_info["longitude"] - self.timezone = site_info["timezone"] - self.elevation = site_info["elevation"] - self.orient = site_info["orientation"] - self.dims = cfg["Dimensions"] - self.filestruct = cfg["FileStructure"] - self.matricesd = os.path.join(self.root, self.filestruct["matrices"]) - self.bsdfd = os.path.join(self.root, self.filestruct["resources"], "BSDFs") - simctrl = cfg["SimulationControl"] - self.parallel = simctrl["parallel"] == "True" - self.vmx_opt = simctrl["vmx"] - self.vsmx_opt = simctrl["vsmx"] - self.dmx_opt = simctrl["dmx"] - self.remake_matrices = simctrl["re-make_matrices"] == "True" - self.view = cfg["View"]["view1"] - self.grid_height = float(cfg["Grid"]["height"]) - self.grid_spacing = float(cfg["Grid"]["spacing"]) - - def get_paths(self): - """Function to collect file paths.""" - - btdfd_vis = os.path.join(self.bsdfd, "vis") - btdfd_shgc = os.path.join(self.bsdfd, "shgc") - btdfd_sol = os.path.join(self.bsdfd, "sol") - self.btdfs_vis = sorted( - [ - os.path.join(self.root, btdfd_vis, i) - for i in os.listdir(btdfd_vis) - if i.endswith(".mtx") - ] - ) - self.btdfs_shgc = sorted( - [ - os.path.join(self.root, btdfd_shgc, i) - for i in os.listdir(btdfd_shgc) - if i.endswith(".mtx") - ] - ) - self.btdfs_tsol = sorted( - [ - os.path.join(self.root, btdfd_sol, i) - for i in os.listdir(btdfd_sol) - if "_Tsol" in i - ] - ) - self.btdfs_ftsol = sorted( - [ - os.path.join(self.root, btdfd_sol, i) - for i in os.listdir(btdfd_sol) - if "fTsol" in i - ] - ) - self.abs1 = sorted( - [ - os.path.join(self.root, btdfd_sol, i) - for i in os.listdir(btdfd_sol) - if "Abs1" in i - ] - ) - self.abs2 = sorted( - [ - os.path.join(self.root, btdfd_sol, i) - for i in os.listdir(btdfd_sol) - if "Abs2" in i - ] - ) - - def make_room(self): - """Function to make a side-lit shoebox room.""" - - windows = [ - [float(n) for n in v.split()] for k, v in self.dims.items() if "window" in k - ] - self.theroom = room.create_south_facing_room( - width=float(self.dims["width"]), - depth=float(self.dims["depth"]), - floor_floor=float(self.dims["height"]) * 1.3, - floor_ceiling=float(self.dims["height"]), - swall_thickness=float(self.dims["facade_thickness"]), - wpd=windows, - ) - self.sensor_grid = utils.gen_grid( - self.theroom.floor.base.flip(), self.grid_height, self.grid_spacing - ) - - def get_workflow(self): - """Function to create workflow.""" - _view = self.view.split() - model = self.theroom.model_dump() - with open("test.rad", "wb") as wtr: - wtr.write(model["materials"]["bytes"]) - wtr.write(model["scene"]["bytes"]) - for w in model["windows"].values(): - wtr.write(w["bytes"]) - model["sensors"] = { - "wpi": {"data": self.sensor_grid}, - "view": {"data": [_view]}, - } - model["views"] = {} - settings = { - "method": 3, - "latitude": self.lat, - "longitude": self.lon, - "time_zone": self.timezone, - "site_elevation": self.elevation, - "save_matrices": True, - "matrix_dir": self.matricesd, - "sensor_window_matrix": self.vmx_opt.split(), - "daylight_matrix": self.dmx_opt.split(), - "sky_basis": "r4", - } - workflow_config = methods.WorkflowConfig.from_dict( - {"model": model, "settings": settings} - ) - self.workflow = methods.ThreePhaseMethod(workflow_config) - - def make_matrices(self): - """Function to make matrices.""" - - self.workflow.generate_matrices(view_matrices=False) - - def load_mtx(self): - """Function to load mtx.""" - self._tvis = [matrix.load_matrix(i) for i in self.btdfs_vis] - self._tsol = [matrix.load_matrix(i) for i in self.btdfs_tsol] - self._ftsol = [matrix.load_matrix(i) for i in self.btdfs_ftsol] - self._abs1 = [matrix.load_matrix(i) for i in self.abs1] - self._abs2 = [matrix.load_matrix(i) for i in self.abs2] - self._shg = [matrix.load_matrix(i) for i in self.btdfs_shgc] - - def wpi_map(self, data, center_w, center_d, fov): - """Function to make wpi map.""" - - height = float(self.dims["height"]) - a = fov / 2 - view = np.tan(np.deg2rad(a)) * (height - self.grid_height) - map_center = data.copy(deep=True) - map_center.index = map_center.index - center_w - map_center.columns = map_center.columns - center_d - map_center = map_center.apply(lambda x: np.sqrt(x.index**2 + x.name**2)) - map_center[map_center > view] = -1 - map_center[(map_center <= view) & (map_center >= 0)] = 0 - map_center = map_center + 1 - map_center = map_center.replace(0, np.nan) - return map_center.values - - def process_wpi(self, wpi_grid): - """Function to process wpi.""" - - if self.init: - width = float(self.dims["width"]) - depth = float(self.dims["depth"]) - - spacing = self.grid_spacing - n_w = int(width / spacing) - n_d = int(depth / spacing) - dist_w = round((width - (n_w - 1) * spacing) / 2, 4) - dist_d = round((depth - (n_d - 1) * spacing) / 2, 4) - - data = pd.DataFrame(wpi_grid.reshape(n_w, n_d)) - self.wpi_data = data.copy(deep=True) - self.wpi_data.index = ( - pd.Series([dist_w] + [spacing for _ in data.index[1:]]).cumsum().values - ) - self.wpi_data.columns = ( - pd.Series([dist_d] + [spacing for _ in data.columns[1:]]) - .cumsum() - .values - ) - - # Center - center_w = round(width / 2, 4) - center_d = round(depth / 2, 4) - fov = 60 # deg - self.map_center = self.wpi_map(self.wpi_data, center_w, center_d, fov) - - # Occupant - occ_w = round(width / 2, 4) - # occ_d = 1.525 # Occupant - occ_d = round(depth / 3 * 2, 4) - fov = 60 # deg - self.map_23 = self.wpi_map(self.wpi_data, occ_w, occ_d, fov) - - self.n_w = n_w - self.n_d = n_d - - data = pd.DataFrame(wpi_grid.reshape(self.n_w, self.n_d)) - data_center = data * self.map_center - wpi_23 = data * self.map_23 - - if self.wpi_plot: - fig, axs = plt.subplots(3, 1, figsize=(6, 10), sharex=False, sharey=False) - tt_avg = round(data.mean().mean(), 1) - axs[0].set_title(f"All WPI sensors (avg={tt_avg})") - pos = axs[0].imshow(data.values, cmap="hot") - axs[0].set_xticks(np.arange(0, self.n_d, 1)) - axs[0].set_xticklabels( - [round(self.wpi_data.columns[int(i)], 1) for i in axs[0].get_xticks()] - ) - axs[0].set_yticks(np.arange(0, self.n_w, 1)) - axs[0].set_yticklabels( - [round(self.wpi_data.index[int(i)], 1) for i in axs[0].get_yticks()] - ) - fig.colorbar(pos, ax=axs[0]) - - tt_avg = round(data_center.mean().mean(), 1) - axs[1].set_title(f"Cone view 60 deg room centered (avg={tt_avg})") - pos = axs[1].imshow(data_center.values, cmap="hot") - axs[1].set_xticks(np.arange(0, self.n_d, 1)) - axs[1].set_xticklabels( - [round(self.wpi_data.columns[int(i)], 1) for i in axs[1].get_xticks()] - ) - axs[1].set_yticks(np.arange(0, self.n_w, 1)) - axs[1].set_yticklabels( - [round(self.wpi_data.index[int(i)], 1) for i in axs[1].get_yticks()] - ) - fig.colorbar(pos, ax=axs[1]) - - tt_avg = round(wpi_23.mean().mean(), 1) - axs[2].set_title(f"Cone view 60 deg occupant centered (avg={tt_avg})") - pos = axs[2].imshow(wpi_23.values, cmap="hot") - axs[2].set_xticks(np.arange(0, self.n_d, 1)) - axs[2].set_xticklabels( - [round(self.wpi_data.columns[int(i)], 1) for i in axs[2].get_xticks()] - ) - axs[2].set_yticks(np.arange(0, self.n_w, 1)) - axs[2].set_yticklabels( - [round(self.wpi_data.index[int(i)], 1) for i in axs[2].get_yticks()] - ) - fig.colorbar(pos, ax=axs[2]) - - fig.tight_layout() - tt_dt = self.date_time.strftime("%Y%m%dT%H%M") - fig.savefig(f"wpi-plot_{tt_dt}.jpeg") - - wpi_all = data.mean().mean() - wpi_center = data_center.mean().mean() - wpi_23 = wpi_23.mean().mean() - - if self.wpi_loc == "23back": - wpi = wpi_23 - elif self.wpi_loc == "center": - wpi = wpi_center - elif self.wpi_loc == "all": - wpi = wpi_all - else: - print(f'ERROR: wpi_loc "{self.wpi_loc}" not available.') - - return wpi, wpi_all, wpi_center, wpi_23 - - def compute(self, weather_df): - """Compute the resulting illuminance and solar gain for all shade states. - Arguments: - weather_df: datetime indexed pandas dataframe with DNI and DHI data - Return: - wpi_df, shg_df, ev_df: Datetime indexed dataframes for each window zone and - window state. wz0-3: bottom-to-top; ws0-n: sorted window state file name from - small to large. - """ - - weather_df = weather_df.copy() - col = weather_df.columns - weather_df.loc[:, "month"] = weather_df.index.month - weather_df.loc[:, "day"] = weather_df.index.day - weather_df.loc[:, "hours"] = ( - weather_df.index.hour + weather_df.index.minute / 60 - ) - _col = ["month", "day", "hours"] - _col.extend(col) - _df = weather_df[_col] - _df = _df.round({"hours": 3}) - sky_data = [" ".join(i.split()[2:]) for i in _df.to_string().splitlines()[1:]] - ntime = len(sky_data) - vissky = makesky.gendaymtx( - sky_data, - self.lat, - self.lon, - self.timezone, - self.elevation, - rotate=self.orient, - )[0] - solsky = makesky.gendaymtx( - sky_data, - self.lat, - self.lon, - self.timezone, - self.elevation, - solar=True, - rotate=self.orient, - )[0] - viscvt = "rmtxop -fa -c 47.4 119.9 11.6" - # solcvt = 'rmtxop -fa -c .265 .67 .065' - shgccvt = "rmtxop -fa -c 1 0 0" - wrapsky = '"' if os.name == "nt" else "'" - cmds = [] - for wname in self.theroom.swall.windows: - wndw_zone = self.theroom.swall.windows[wname] - for idx_t, btdf_vis in enumerate(self.btdfs_vis): - vis_cmd = f"rmtxop {self.vmxs[wname]} {btdf_vis} {self.dmxs[wname]}" - vis_cmd += f" {wrapsky}!{vissky}{wrapsky} | {viscvt} - | getinfo - " - shg_cmd = f"rmtxop {self.btdfs_shgc[idx_t]} {self.dmxs[wname]}" - shg_cmd += f" {wrapsky}!{solsky}{wrapsky} | {shgccvt} - |" - shg_cmd += f" rmtxop -fa -s {wndw_zone.area()} - | getinfo -" - if os.name == "nt": - vis_cmd = self.skycmd_win(vis_cmd) - shg_cmd = self.skycmd_win(shg_cmd) - cmds.append(vis_cmd) - cmds.append(shg_cmd) - if self.parallel: - with mp.Pool(mp.cpu_count()) as pool: - raw = [i.decode().split() for i in pool.map(checkout, cmds)] - else: - # raw = [checkout(cmd).decode().split() for cmd in cmds] - - raw = [] - for cmd in cmds: - st = time.time() - raw.append(checkout(cmd).decode().split()) - et = time.time() - print(cmd[:10], round(et - st, 2)) - - col = [ - f"{zone}_{state}" - for zone in range(len(self.theroom.swall.windows)) - for state in range(len(self.btdfs_vis)) - ] - vis_df = pd.DataFrame(raw[::2], dtype=float).transpose() - wpi_df = vis_df.iloc[:-ntime] - awpi_df = pd.DataFrame( - [wpi_df.iloc[n::ntime].mean() for n in range(ntime)], index=_df.index - ) - awpi_df.columns = ["wpi_" + c for c in col] - ev_df = vis_df.iloc[-ntime:].set_index(_df.index) - ev_df.columns = ["vil_" + c for c in col] - shg_df = pd.DataFrame(raw[1::2], dtype=float).transpose().set_index(_df.index) - shg_df.columns = ["shg_" + c for c in col] - return pd.concat([awpi_df, shg_df, ev_df], axis=1) - - def compute2(self, weather_df, weather_cutoff=75): - """Compute the resulting illuminance and solar gain for all shade states. - Arguments: - weather_df: datetime indexed pandas dataframe with DNI and DHI data - weather_cutoff: cutoff for weather forecast, W/m2 - Return: - output_df: Datetime indexed dataframes for each window zone and - window state. wz0-3: bottom-to-top; ws0-n: sorted window state file name from - small to large. - """ - - weather_df = weather_df.copy() - - # Filter weather - wather_filter = (weather_df["DNI"] <= weather_cutoff) & ( - weather_df["DHI"] <= weather_cutoff - ) - weather_df["DNI"] = weather_df["DNI"].mask(wather_filter, 0) - weather_df["DHI"] = weather_df["DHI"].mask(wather_filter, 0) - - col = weather_df.columns - weather_df.loc[:, "month"] = weather_df.index.month - weather_df.loc[:, "day"] = weather_df.index.day - weather_df.loc[:, "hours"] = ( - weather_df.index.hour + weather_df.index.minute / 60 - ) - _col = ["month", "day", "hours"] - _col.extend(col) - _df = weather_df[_col] - _df = _df.round({"hours": 3}) - - _vis = self.workflow.get_sky_matrix( - time=list(weather_df.index.tolist()), - dni=weather_df["DNI"].tolist(), - dhi=weather_df["DHI"].tolist(), - ) - _sol = self.workflow.get_sky_matrix( - time=list(weather_df.index.tolist()), - dni=weather_df["DNI"].tolist(), - dhi=weather_df["DHI"].tolist(), - solar_spectrum=True, - ) - output_df = pd.DataFrame(index=_df.index) - flr_area = np.linalg.norm(self.theroom.floor.base.area) - for idx, window in enumerate(self.theroom.swall.windows): - _wndw_area = np.linalg.norm(window.polygon.area) - _window_name = window.primitive.identifier - for idx_t, _ in enumerate(self.btdfs_vis): - wpi_col = f"wpi_{idx}_{idx_t}" - ev_col = f"ev_{idx}_{idx_t}" - shg_col = f"shg_{idx}_{idx_t}" - tsol_col = f"tsol_{idx}_{idx_t}" - abs1_col = f"abs1_{idx}_{idx_t}" - abs2_col = f"abs2_{idx}_{idx_t}" - iflr_col = f"iflr_{idx}_{idx_t}" - visres = matrix.matrix_multiply_rgb( - self.workflow.sensor_window_matrices["wpi"].array[idx], - self._tvis[idx_t], - self.workflow.daylight_matrices[_window_name].array, - _vis, - weights=[47.4, 119.9, 11.6], - ) - evres = matrix.matrix_multiply_rgb( - self.workflow.sensor_window_matrices["view"].array[0], - self._tvis[idx_t], - self.workflow.daylight_matrices[_window_name].array, - _vis, - weights=[47.4, 119.9, 11.6], - ) - - # Process WPI - wpi = [] - for i, t in enumerate(np.transpose(visres)): - self.date_time = weather_df.index[i] - wpi_temp = self.process_wpi(t) - wpi.append(wpi_temp[0]) - output_df[wpi_col] = wpi - output_df[ev_col] = evres[0] - shg = matrix.matrix_multiply_rgb( - self._shg[idx_t], - self.workflow.daylight_matrices[_window_name].array, - _sol, - weights=[.3, .4, .3], - )[0] - abs1 = matrix.matrix_multiply_rgb( - self._abs1[idx_t], - self.workflow.daylight_matrices[_window_name].array, - _sol, - weights=[.3, .4, .3], - )[0] - abs2 = matrix.matrix_multiply_rgb( - self._abs2[idx_t], - self.workflow.daylight_matrices[_window_name].array, - _sol, - weights=[.3, .4, .3], - )[0] - tsol = matrix.matrix_multiply_rgb( - self._ftsol[idx_t], - self.workflow.daylight_matrices[_window_name].array, - _sol, - weights=[.3, .4, .3], - ) - iflr = matrix.matrix_multiply_rgb( - self.workflow.surface_window_matrices["floor"].array[idx], - self._ftsol[idx_t], - self.workflow.daylight_matrices[_window_name].array, - _sol, - weights=[.3, .4, .3], - )[0] - output_df[shg_col] = shg * _wndw_area - output_df[tsol_col] = (tsol.T * KFOMG).sum(axis=1) * _wndw_area - output_df[iflr_col] = iflr * np.pi * flr_area - output_df[abs1_col] = abs1 * _wndw_area - output_df[abs2_col] = abs2 * _wndw_area - - # Process if shade - if "shade" in self.facade_type and not self._test_difference: - if not self.new_map: - cols = output_df.columns - for prefix in np.unique([c.split("_")[0] for c in cols]): - for i, v in enumerate(list(self.shade_map.values())[:-1]): - self.new_map[f"{prefix}_0_{i}"] = [ - f"{prefix}_{ii}_{vv}" for ii, vv in enumerate(v) - ] - i += 1 - self.new_map[f"{prefix}_0_{i}"] = [ - f"{prefix}_{ii}_1" for ii, vv in enumerate(v) - ] - temp = output_df.copy(deep=True) - for k, v in self.new_map.items(): - output_df[k] = temp[v].sum(axis=1) - output_df = output_df[self.new_map.keys()] - - self.init = False - - return output_df - - -if __name__ == "__main__": - - import warnings - warnings.simplefilter(action='ignore', category=pd.errors.PerformanceWarning) - - from .configs import get_config - root = os.path.dirname(os.path.abspath(__file__)) - - # configuration - wwr = 0.6 # [0.4, 0.6] - mode = 'dshade' # ['shade', 'blinds', 'ec'] - facade_type = 'shade' # ['shade', 'blinds', 'ec'] - print('Running example for:', wwr, mode, facade_type) - filestruct, config_path = get_config(mode, wwr) - forecaster = Forecast( - config_path, - facade_type=facade_type, - regenerate=False, - filestruct=filestruct, - wpi_plot=False, - wpi_loc="23back", - ) - - for i in range(3): - #df = pd.DataFrame([[800, 80], [750, 100], [100, 10]], - # index=pd.DatetimeIndex([ - # pd.datetime(2019, 12, 21, 12, 0), - # pd.datetime(2019, 12, 21, 12, 10), - # pd.datetime(2019, 12, 21, 12, 20), - # ])) - df = pd.DataFrame(index=pd.date_range('2020-01-01 00:00', - '2020-01-02 00:00', - freq='5T')) - np.random.seed(1) - df['DNI'] = np.random.uniform(0, 1000, len(df)) - np.random.seed(1) - df['DHI'] = np.random.uniform(0, 250, len(df)) - - fake_df = pd.DataFrame([[185.5, 42], [165.3, 38.3]], - index=pd.DatetimeIndex([pd.to_datetime('2017-09-01 18:15'), - pd.to_datetime('2017-09-01 18:20')]), - columns=['DNI','DHI']) - - st = time.time() - #res = forecaster.compute(df) - res = forecaster.compute2(fake_df) - print(time.time() - st) - print(f'Forecast of 24h x 5min in {round(time.time() - st, 2)} s') - #print(res.columns) - #for k,v in forecaster.new_map.items(): - # print(k, v) - print(res) - print(res['wpi_0_6']) - res.to_csv('radiance-forecast_new.csv') diff --git a/afc/resources/radiance/matrices/shade/0.6/d654a92a0aef00f4.npz b/afc/resources/radiance/matrices/shade/0.6/d654a92a0aef00f4.npz deleted file mode 100644 index 5ba0c1e..0000000 Binary files a/afc/resources/radiance/matrices/shade/0.6/d654a92a0aef00f4.npz and /dev/null differ