28 11

Author: paulaharder

evaluate.py

@@ -51,6 +51,7 @@ def add_arguments():
  parser.add_argument("--factor", type=int, default=2)
  parser.add_argument("--time_sr", default=False)
  parser.add_argument("--constraints_window_size", default=4, type=int)
+ parser.add_argument("--ensemble", default=False)
 
  return parser.parse_args()
 

mean_scores.py

@@ -0,0 +1,43 @@
+import csv
+
+def add_arguments():
+ parser = argparse.ArgumentParser()
+ parser.add_argument("--model_id", default='test')
+ return parser.parse_args()
+
+#load five csvs and calculate means and stds
+
+def main(args):
+ number_runs = 5
+ metrics = ['MSE', 'RMSE', 'PSNR', 'MAE', 'SSIM', 'Mass_violation', 'Mean bias', 'MS SSIM', 'Pearson corr', 'CRPS']
+ dict_lists = {}
+ for metric in metrics:
+ means[metric] = []
+ for i in range(number_runs):#
+ filename = args.model_id + '_' + str(i)
+ with open(filename,'r') as data:
+ for line in csv.reader(data):
+ if line[0] in metrics:
+ means[line[0]].append(line[1])
+ 
+ #iterate over dict lists
+ for metric, values in state_dict.items():
+ dict_lists[metric+'_mean'] = np.mean(np.array(values))
+ dict_lists[metric+'_std'] = np.std(np.array(values))
+ 
+ #save mean+std dict as csv
+ save_dict(dict_lists, args)
+
+ 
+ 
+def save_dict(dictionary, args):
+ w = csv.writer(open('./data/score_log/'+args.model_id+'_means.csv', 'w'))
+ # loop over dictionary keys and values
+ for key, val in dictionary.items():
+ # write every key and value to file
+ w.writerow([key, val])
+ 
+
+if __name__ == '__main__':
+ args = add_arguments()
+ main(args)

models.py

@@ -434,6 +434,7 @@ def __init__(self, upsampling_factor):
  self.pool = torch.nn.AvgPool2d(kernel_size=upsampling_factor)
  self.upsampling_factor = upsampling_factor
  def forward(self, y, lr):
+ y = y.clone()
  out = self.pool(y)
  out = y*torch.kron(lr*1/out, torch.ones((self.upsampling_factor,self.upsampling_factor)).to('cuda'))
  return out
@@ -700,10 +701,84 @@ def forward(self, x, mr=None, z=None):
  # out[:,0,i,...] = self.constraints(out, x[:,0,i,...])
  #out[:,0,:,:] *= 16
  out = out.unsqueeze(1)
+ return out
+ 
+class ResNet3(nn.Module):
+ def __init__(self, number_channels=64, number_residual_blocks=4, upsampling_factor=2, noise=False, constraints='none', dim=1, cwindow_size=2):
+ super(ResNet3, self).__init__()
+ # First layer
+ if noise:
+ self.conv_trans0 = nn.ConvTranspose2d(100, 1, kernel_size=(32,32), padding=0, stride=1)
+ self.conv1 = nn.Sequential(nn.Conv2d(dim, number_channels, kernel_size=3, stride=1, padding=1), nn.ReLU(inplace=True))
+ else:
+ self.conv1 = nn.Sequential(nn.Conv2d(dim, number_channels, kernel_size=3, stride=1, padding=1), nn.ReLU(inplace=True))
+ #Residual Blocks
+ self.res_blocks = nn.ModuleList()
+ for k in range(number_residual_blocks):
+ self.res_blocks.append(ResidualBlock(number_channels, number_channels))
+ # Second conv layer post residual blocks
+ self.conv2 = nn.Sequential(
+ nn.Conv2d(number_channels, number_channels, kernel_size=3, stride=1, padding=1), nn.ReLU(inplace=True))
+ # Upsampling layers
+ self.upsampling = nn.ModuleList()
+ for k in range(1):
+ self.upsampling.append(nn.ConvTranspose2d(number_channels, number_channels, kernel_size=3, padding=0, stride=3) )
+ # Next layer after upper sampling
+ self.conv3 = nn.Sequential(nn.Conv2d(number_channels, number_channels, kernel_size=3, stride=1, padding=1), nn.ReLU(inplace=True))
+ # Final output layer
+ self.conv4 = nn.Conv2d(number_channels, dim, kernel_size=1, stride=1, padding=0) 
+ #optional renomralization layer
+ self.is_constraints = False
+ if constraints == 'softmax':
+ self.constraints = SoftmaxConstraints(upsampling_factor=upsampling_factor, cwindow_size=cwindow_size)
+ self.is_constraints = True
+ elif constraints == 'enforce_op':
+ self.constraints = EnforcementOperator(upsampling_factor=upsampling_factor)
+ self.is_constraints = True
+ elif constraints == 'add':
+ self.constraints = AddDownscaleConstraints(upsampling_factor=upsampling_factor)
+ self.is_constraints = True
+ elif constraints == 'mult':
+ self.constraints = MultDownscaleConstraints(upsampling_factor=upsampling_factor)
+ self.is_constraints = True
+ 
+ self.dim = dim 
+ self.noise = noise
+ 
+ def forward(self, x, mr=None, z=None): 
+ if self.noise:
+ out = self.conv_trans0(z)
+ out = self.conv1(torch.cat(( x[:,0,...],out), dim=1))
+ for layer in self.res_blocks:
+ out = layer(out)
+ out = self.conv2(out)
+ for layer in self.upsampling:
+ out = layer(out)
+ out = self.conv3(out)
+ out = self.conv4(out)
+ if self.is_constraints:
+ out = self.constraints(out, x[:,0,...])
  return out 
+ else:
+ #print(x.shape)
+ out = self.conv1(x[:,0,...])
+ for layer in self.upsampling:
+ out = layer(out)
+ out = self.conv2(out) 
+ for layer in self.res_blocks:
+ out = layer(out)
+ out = self.conv3(out)
+ out = self.conv4(out)
+ if self.is_constraints:
+ out[:,...] = self.constraints(out, x[:,0,...])
+ #for i in range(self.dim):
+ # out[:,0,i,...] = self.constraints(out, x[:,0,i,...])
+ #out[:,0,:,:] *= 16
+ out = out.unsqueeze(1)
+ return out
 
 class ResNetNoise(nn.Module):
- def __init__(self, number_channels=64, number_residual_blocks=4, upsampling_factor=2, noise=False, constraints='none', dim=1):
+ def __init__(self, number_channels=64, number_residual_blocks=4, upsampling_factor=2, noise=False, constraints='none', dim=1, cwindow_size=2):
  super(ResNetNoise, self).__init__()
  # First layer
 
@@ -728,7 +803,7 @@ def __init__(self, number_channels=64, number_residual_blocks=4, upsampling_fact
  #optional renomralization layer
  self.is_constraints = False
  if constraints == 'softmax':
- self.constraints = SoftmaxConstraints(upsampling_factor=upsampling_factor)
+ self.constraints = SoftmaxConstraints(upsampling_factor=upsampling_factor, cwindow_size=cwindow_size)
  self.is_constraints = True
  elif constraints == 'enforce_op':
  self.constraints = EnforcementOperator(upsampling_factor=upsampling_factor)
@@ -817,7 +892,7 @@ def forward(self, x):
 
 
 class ResNet2Up(nn.Module):
- def __init__(self, number_channels=64, number_residual_blocks=4, upsampling_factor=2, noise=False, constraints='none', dim=1, output_mr=False):
+ def __init__(self, number_channels=64, number_residual_blocks=4, upsampling_factor=2, noise=False, constraints='none', dim=1, output_mr=False, cwindow_size= 2):
  super(ResNet2Up, self).__init__()
  #PART I
  self.conv1 = nn.Sequential(nn.Conv2d(dim, number_channels, kernel_size=3, stride=1, padding=1), nn.ReLU(inplace=True))
@@ -833,7 +908,7 @@ def __init__(self, number_channels=64, number_residual_blocks=4, upsampling_fact
 
  self.is_constraints = False
  if constraints == 'softmax':
- self.constraints = SoftmaxConstraints(upsampling_factor=2)
+ self.constraints = SoftmaxConstraints(upsampling_factor=2, cwindow_size=cwindow_size)
  self.is_constraints = True
  elif constraints == 'enforce_op':
  self.constraints = EnforcementOperator(upsampling_factor=2)
@@ -844,7 +919,7 @@ def __init__(self, number_channels=64, number_residual_blocks=4, upsampling_fact
 
  self.is_single_constraints = False
  if constraints == 'softmax_single':
- self.constraints = SoftmaxConstraints(upsampling_factor=4)
+ self.constraints = SoftmaxConstraints(upsampling_factor=4, cwindow_size=cwindow_size)
  self.is_single_constraints = True
  elif constraints == 'enforce_op_single':
  self.constraints = EnforcementOperator(upsampling_factor=4)
@@ -865,7 +940,7 @@ def __init__(self, number_channels=64, number_residual_blocks=4, upsampling_fact
  self.conv23 = nn.Sequential(nn.Conv2d(number_channels, number_channels, kernel_size=3, stride=1, padding=1), nn.ReLU(inplace=True))
  self.conv24 = nn.Conv2d(number_channels, dim, kernel_size=1, stride=1, padding=0)
  if constraints == 'softmax':
- self.constraints2 = SoftmaxConstraints(upsampling_factor=2)
+ self.constraints2 = SoftmaxConstraints(upsampling_factor=2, cwindow_size=cwindow_size)
  elif constraints == 'enforce_op':
  self.constraints2 = EnforcementOperator(upsampling_factor=2)
  elif constraints == 'mult':

scoring.py

@@ -7,6 +7,8 @@
 import argparse
 from torch.utils.data import DataLoader, TensorDataset
 import csv
+from torchmetrics.functional import multiscale_structural_similarity_index_measure, structural_similarity_index_measure
+
 
 def add_arguments():
  parser = argparse.ArgumentParser()
@@ -20,28 +22,40 @@ def add_arguments():
  parser.add_argument("--mass_violation", type=bool, default=True)
  parser.add_argument("--factor", type=int, default=4)
  parser.add_argument("--time_sr", default=False)
+ 
  #args for model loading
  return parser.parse_args()
 
 def main_scoring(args):
  #n = 24
  input_val = torch.load('./data/test/'+args.dataset+'/input_test.pt')
  target_val = torch.load('./data/test/'+args.dataset+'/target_test.pt')
+ #target_val = torch.load('./data/test/'+args.dataset+'/target_test.pt')
  #target_val = torch.load('./data/test/dataset28/target_test.pt')
  val_data = TensorDataset(input_val, target_val)
  pred = np.zeros(target_val.shape)
+ max_val = target_val.max()
+ min_val = target_val.min()
  print(pred.shape)
  factor = args.factor
  mse = 0
  mae = 0
  ssim = 0
+ mean_bias = 0
  mass_violation = 0
+ ms_ssim = 0
+ corr = 0
+ crps = 0
+ 
  l2_crit = nn.MSELoss()
  l1_crit = nn.L1Loss()
- ssim_criterion = tgm.losses.SSIM(window_size=11, max_val=130.83, reduction='mean')
+ #ssim_criterion = StructuralSimilarityIndexMeasure() #tgm.losses.SSIM(window_size=11, max_val=max_val, reduction='mean')
  if args.nn:
- 
- pred = torch.load('./data/prediction/'+args.dataset+'_'+args.model_id+'_test.pt')
+ if args.ensemble:
+ en_pred = torch.load('./data/prediction/'+args.dataset+'_'+args.model_id+'_test_ensemble.pt')
+ pred = torch.mean(en_pred, dim=1)
+ else:
+ pred = torch.load('./data/prediction/'+args.dataset+'_'+args.model_id+'_test.pt')
  pred = pred.detach().cpu().numpy()
  print(pred.shape)
  for i,(lr, hr) in enumerate(val_data):
@@ -71,7 +85,13 @@ def main_scoring(args):
  print(i, mse_loss)'''
  mse += l2_crit(torch.Tensor(pred[i,j,...]), hr[j,...]).item()
  mae += l1_crit(torch.Tensor(pred[i,j,...]), hr[j,...]).item()
- ssim += ssim_criterion(torch.Tensor(pred[i,j:j+1,...]), hr[j:j+1,...]).item()
+ mean_bias += torch.mean( hr[j,...]-torch.abs(torch.Tensor(pred[i,j,...])))
+ corr += pearsonr(torch.Tensor(pred[i,j,...]).flatten(), hr[j,...].flatten())
+ ms_ssim += multiscale_structural_similarity_index_measure(torch.Tensor(pred[i,j:j+1,...]), hr[j:j+1,...], data_range=max_val-min_val, kernel_size=11, betas=(0.2856, 0.3001, 0.2363))#0.0448, 0.2856, 0.3001, 0.2363))
+ ssim += structural_similarity_index_measure(torch.Tensor(pred[i,j:j+1,...]), hr[j:j+1,...] , data_range=max_val-min_val, kernel_size=11)#ssim_criterion(torch.Tensor(pred[i,j:j+1,...]), hr[j:j+1,...]).item()
+ if args.ensemble:
+ crps_ens = crps_ensemble(hr[j,0,0,...].numpy(), np.swapaxis(np.swapaxis(pred[i,:,j,0,0,...], 0,1),1,2))
+ crps += np.mean(crps_ens)
  if args.mass_violation:
  if args.time_sr:
  if j==0:
@@ -105,24 +125,81 @@ def main_scoring(args):
  mse *= 1/(input_val.shape[0]*args.time_steps)
  mae *= 1/(input_val.shape[0]*args.time_steps)
  ssim *= 1/(input_val.shape[0]*args.time_steps)
+ mean_bias *= 1/(input_val.shape[0]*args.time_steps)
+ corr *= 1/(input_val.shape[0]*args.time_steps)
+ ms_ssim *= 1/(input_val.shape[0]*args.time_steps)
+ 
  if args.mass_violation:
  if args.time_sr:
  mass_violation *= 1/(input_val.shape[0]*args.time_steps)
  else:
- mass_violation *= 1/(input_val.shape[0]*2)
+ mass_violation *= 1/(input_val.shape[0]) #what is the 2 doing here?
  else:
  mse *= 1/input_val.shape[0] 
  mae *= 1/input_val.shape[0] 
  ssim *= 1/input_val.shape[0] 
  psnr = calculate_pnsr(mse, target_val.max() ) 
- scores = {'MSE':mse, 'RMSE':torch.sqrt(torch.Tensor([mse])), 'PSNR': psnr, 'MAE':mae, 'SSIM':1-ssim, 'Mass_violation': mass_violation}
+ scores = {'MSE':mse, 'RMSE':torch.sqrt(torch.Tensor([mse])), 'PSNR': psnr, 'MAE':mae, 'SSIM':ssim, 'Mass_violation': mass_violation, 'Mean bias': mean_bias, 'MS SSIM': ms_ssim, 'Pearson corr': corr, 'CRPS': crps}
  print(scores)
  create_report(scores, args)
  #np.save('./data/prediction/bic.npy', pred)
 
 
 def calculate_pnsr(mse, max_val):
  return 20 * torch.log10(max_val / torch.sqrt(torch.Tensor([mse])))
+ 
+def pearsonr(x, y):
+ """
+ Mimics `scipy.stats.pearsonr`
+ Arguments
+ ---------
+ x : 1D torch.Tensor
+ y : 1D torch.Tensor
+ Returns
+ -------
+ r_val : float
+ pearsonr correlation coefficient between x and y
+ 
+ Scipy docs ref:
+ https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html
+ 
+ Scipy code ref:
+ https://github.com/scipy/scipy/blob/v0.19.0/scipy/stats/stats.py#L2975-L3033
+ Example:
+ >>> x = np.random.randn(100)
+ >>> y = np.random.randn(100)
+ >>> sp_corr = scipy.stats.pearsonr(x, y)[0]
+ >>> th_corr = pearsonr(torch.from_numpy(x), torch.from_numpy(y))
+ >>> np.allclose(sp_corr, th_corr)
+ """
+ mean_x = torch.mean(x)
+ mean_y = torch.mean(y)
+ xm = x.sub(mean_x)
+ ym = y.sub(mean_y)
+ r_num = xm.dot(ym)
+ r_den = torch.norm(xm, 2) * torch.norm(ym, 2)
+ r_val = r_num / r_den
+ return r_val
+
+def crps_ensemble(observation, forecasts):
+ fc = forecasts.copy()
+ fc.sort(axis=-1)
+ obs = observation
+ fc_below = fc<obs[...,None]
+ crps = np.zeros_like(obs)
+
+ for i in range(fc.shape[-1]):
+ below = fc_below[...,i]
+ weight = ((i+1)**2 - i**2) / fc.shape[-1]**2
+ crps[below] += weight * (obs[below]-fc[...,i][below])
+
+ for i in range(fc.shape[-1]-1,-1,-1):
+ above = ~fc_below[...,i]
+ k = fc.shape[-1]-1-i
+ weight = ((k+1)**2 - k**2) / fc.shape[-1]**2
+ crps[above] += weight * (fc[...,i][above]-obs[above])
+
+ return crps
 
 def create_report(scores, args):
  args_dict = args_to_dict(args)