add

Pk_analysis/.ipynb_checkpoints/analize_study-checkpoint.py

@@ -0,0 +1,71 @@
+# This script load a study and output some information about it
+import numpy as np
+import optuna
+
+################################## INPUT ########################################
+# study parameters
+study_name = 'params_2_SFRH_TPE'
+storage    = 'sqlite:///TPE.db'
+
+# whether print information of a particular trial
+individual_trials = None #if None, nothing is printed
+#################################################################################
+
+# load the study
+study = optuna.load_study(study_name=study_name, storage=storage)
+
+# get the number of trials
+trials = len(study.trials)
+print('Found %d trials'%trials)
+
+# get the scores of the trials and print the top-10
+values = np.zeros(trials)
+for i,t in enumerate(study.trials):
+    values[i] = t.value
+
+indexes = np.argsort(values)
+print('\nBest trials:')
+for i in range(10):
+    print('study %04d ----> score: %.5e'%(indexes[i], values[indexes[i]]))
+
+# get the info of the best trial
+trial = study.best_trial
+print("\nBest trial:  number {}".format(trial.number))
+print("  Value: ", trial.value)
+print(" Params: ")
+for key, value in trial.params.items():
+    print("    {}: {}".format(key, value))
+
+if individual_trials is not None:
+    for trial_number in individual_trials:
+        trial = study.trials[trial_number]
+        print(" \nTrial number {}".format(trial_number))
+        print("  Value: ", trial.value)
+        print(" Params: ")
+        for key, value in trial.params.items():
+            print("    {}: {}".format(key, value))
+
+"""
+# read some parameters from the trials and save results to a .txt file
+lr   = np.zeros(10000, dtype=np.float64)
+wd   = np.zeros(10000, dtype=np.float64)
+nl   = np.zeros(10000, dtype=np.float64)
+loss = np.zeros(10000, dtype=np.float64)
+for i,t in enumerate(study.trials):
+    if i>=10000:  break
+    lr[i]   = t.params['lr']
+    wd[i]   = t.params['wd']
+    nl[i]   = t.params['n_layers']
+    loss[i] = t.value
+np.savetxt('borrar.txt',np.transpose([np.arange(10000),nl,lr,wd,loss]))
+"""
+"""
+# get the importances of the different hyperparameters
+importances = optuna.importance.get_param_importances(study)
+print('\nImportances:')
+for key in importances:
+    print('{} {}'.format(key, importances[key]))
+"""
+
+
+

Pk_analysis/.ipynb_checkpoints/architecture-checkpoint.py

@@ -0,0 +1,37 @@
+import torch 
+import torch.nn as nn
+import numpy as np
+import sys, os, time
+import optuna
+
+# This routine returns an architecture that is built inside the routine itself
+# It can have from 1 to max_layers hidden layers. The user specifies the size of the
+# input and output together with the maximum number of neurons in each layers
+# trial -------------> optuna variable
+# input_size --------> size of the input
+# output_size -------> size of the output
+# max_layers --------> maximum number of hidden layers to consider (default=3)
+# max_neurons_layer -> the maximum number of neurons a layer can have (default=500)
+def dynamic_model(trial, input_size, output_size, max_layers=3, max_neurons_layers=500):
+
+    # define the tuple containing the different layers
+    layers = []
+
+    # get the number of hidden layers
+    n_layers = trial.suggest_int("n_layers", 1, max_layers)
+
+    # get the hidden layers
+    in_features = input_size
+    for i in range(n_layers):
+        out_features = trial.suggest_int("n_units_l{}".format(i), 4, max_neurons_layers)
+        layers.append(nn.Linear(in_features, out_features))
+        layers.append(nn.LeakyReLU(0.2))
+        p = trial.suggest_float("dropout_l{}".format(i), 0.2, 0.8)
+        layers.append(nn.Dropout(p))
+        in_features = out_features
+
+    # get the last layer
+    layers.append(nn.Linear(out_features, output_size))
+
+    # return the model
+    return nn.Sequential(*layers)

Pk_analysis/.ipynb_checkpoints/data-checkpoint.py

@@ -0,0 +1,72 @@
+import torch 
+from torch.utils.data.dataset import Dataset
+from torch.utils.data import DataLoader
+import numpy as np
+import sys, os, time
+
+# This class creates the dataset 
+class make_dataset():
+
+    def __init__(self, mode, seed, f_Pk, f_Pk_norm, f_params):
+
+        # read data, scale it, and normalize it
+        Pk = np.log10(np.load(f_Pk))
+        if f_Pk_norm is None:
+            mean, std = np.mean(Pk, axis=0), np.std(Pk, axis=0)
+        else:
+            Pk_norm = np.log10(np.load(f_Pk_norm))
+            mean, std = np.mean(Pk_norm, axis=0), np.std(Pk_norm, axis=0)
+        Pk = (Pk - mean)/std
+
+        # read the value of the cosmological & astrophysical parameters; normalize them
+        params  = np.loadtxt(f_params)
+        params  = np.repeat(params, 24, axis=0)
+
+        minimum = np.array([0.1, 0.6, 0.5, -1.0, -1.0, 80.])
+        maximum = np.array([0.5, 1.0, 2.0,  1.0,  1.0, 1000.])
+        params  = (params - minimum)/(maximum - minimum)
+
+        # get the size and offset depending on the type of dataset
+        sims = Pk.shape[0]
+        if   mode=='train':  size, offset = int(sims*0.90), int(sims*0.00)
+        elif mode=='valid':  size, offset = int(sims*0.05), int(sims*0.90)
+        elif mode=='test':   size, offset = int(sims*0.05), int(sims*0.95)
+        elif mode=='all':    size, offset = int(sims*1.00), int(sims*0.00)
+        else:                raise Exception('Wrong name!')
+
+        # randomly shuffle the sims. Instead of 0 1 2 3...999 have a 
+        # random permutation. E.g. 5 9 0 29...342
+        np.random.seed(seed)
+        indexes = np.arange(sims) #only shuffle realizations, not rotations
+        np.random.shuffle(indexes)
+        indexes = indexes[offset:offset+size] #select indexes of mode
+
+        # select the data in the considered mode
+        Pk     = Pk[indexes]
+        params = params[indexes]
+
+        # define size, input and output matrices
+        self.size   = size
+        self.input  = torch.tensor(Pk,     dtype=torch.float)
+        self.output = torch.tensor(params, dtype=torch.float)
+
+    def __len__(self):
+        return self.size
+
+    def __getitem__(self, idx):
+        return self.input[idx], self.output[idx]
+
+
+# This routine creates a dataset loader
+# mode ---------------> 'train', 'valid', 'test' or 'all'
+# seed ---------------> random seed to split data among training, validation and testing
+# f_Pk ---------------> file containing the power spectra
+# f_Pk_norm ----------> file containing the power spectra to normalize data
+# f_params -----------> files with the value of the cosmological + astrophysical params
+# batch_size ---------> batch size
+# shuffle ------------> whether to shuffle the data or not
+# workers --------> number of CPUs to load the data in parallel
+def create_dataset(mode, seed, f_Pk, f_Pk_norm, f_params, batch_size, shuffle, workers):
+    data_set = make_dataset(mode, seed, f_Pk, f_Pk_norm, f_params)
+    return DataLoader(dataset=data_set, batch_size=batch_size, shuffle=shuffle,
+                      num_workers=workers)

Pk_analysis/.ipynb_checkpoints/gpus_job_Pk2D_9503587-checkpoint.out

@@ -0,0 +1,13 @@
+
+The following have been reloaded with a version change:
+  1) cuda/11.4.4 => cuda/11.4.4-ldlywt5
+
+is available:  True
+device count:  1
+current device:  0
+cuda device:  <torch.cuda.device object at 0x153ea52d97d0>
+cuda device name:  NVIDIA A100-SXM4-80GB
+cuda version:  11.8
+[I 2024-02-05 16:14:37,156] Using an existing study with name 'Pk2D_noise_kmax02_all_params' instead of creating a new one.
+/home/aobulj/data/conda/envs/torch_env/lib/python3.11/site-packages/torch/utils/data/dataloader.py:557: UserWarning: This DataLoader will create 10 worker processes in total. Our suggested max number of worker in current system is 1, which is smaller than what this DataLoader is going to create. Please be aware that excessive worker creation might get DataLoader running slow or even freeze, lower the worker number to avoid potential slowness/freeze if necessary.
+  warnings.warn(_create_warning_msg(

Pk_analysis/.ipynb_checkpoints/logfile_9503586-checkpoint.log

@@ -0,0 +1 @@
+CUDA Available

Pk_analysis/.ipynb_checkpoints/script-checkpoint.sh

@@ -0,0 +1,15 @@
+#!/usr/bin/env bash
+#SBATCH --gpus=1
+#SBATCH --mem=10G
+#SBATCH --constraint="GPUMEM32GB|GPUMEM80GB"
+#SBATCH --time=24:00:00
+#SBATCH [email protected]
+#SBATCH --mail-type=ALL 
+#SBATCH --output=gpus_job_Pk2D_%j.out
+
+module load multigpu cuda/11.4.4
+module load multigpu cudnn/8.2.4
+module load mamba
+source activate torch_env
+python -c 'import torch as t; print("is available: ", t.cuda.is_available()); print("device count: ", t.cuda.device_count()); print("current device: ", t.cuda.current_device()); print("cuda device: ", t.cuda.device(0)); print("cuda device name: ", t.cuda.get_device_name(0)); print("cuda version: ", t.version.cuda)'
+python main.py >> logfile_${SLURM_JOB_ID}.log

Pk_analysis/.ipynb_checkpoints/test-checkpoint.py

@@ -0,0 +1,128 @@
+import numpy as np
+import sys, os, time
+import torch
+import torch.nn as nn
+import data, architecture
+import optuna
+
+################################### INPUT ############################################
+# data parameters
+path        = '/shares/stadel.ics.mnf.uzh/aobulj/hifi_mocks_Paco_noise/' #output folder
+nsim      = 1000
+kmax      = 0.2
+f_Pk      = path + 'Pk2D_maps_%i_kmax2D_%.2f.npy'%(nsim, kmax)
+f_Pk_norm = None
+f_params  = '/home/aobulj/data/Hi-Fi_mocks_Paco/analyse/sobol_sequence_noise_1000.txt'
+seed      = 1                                    #seed to split data in train/valid/test
+mode      = 'test'     #'train','valid','test' or 'all'
+
+# architecture parameters
+input_size  = 31 #number of bins in Pk
+output_size = 12  #number of parameters to predict (posterior mean + std)
+
+# training parameters
+batch_size = 32
+workers    = 1
+g          = [0,1,2,3,4,5]  #minimize loss using parameters
+h          = [6,7,8,9,10,11]  #minimize loss using errors of parameters
+
+# optuna parameters
+study_name = 'Pk2D_noise_kmax02_all_params'
+storage    = 'sqlite:///databases_%s.db'%study_name
+
+# output file name
+fout = 'Results.txt'
+######################################################################################
+
+# use GPUs if available
+if torch.cuda.is_available():
+    print("CUDA Available")
+    device = torch.device('cuda')
+else:
+    print('CUDA Not Available')
+    device = torch.device('cpu')
+
+# load the optuna study
+study = optuna.load_study(study_name=study_name, storage=storage)
+
+# get the scores of the study trials
+values = np.zeros(len(study.trials))
+completed = 0
+for i,t in enumerate(study.trials):
+    values[i] = t.value
+    if t.value is not None:  completed += 1
+
+# get the info of the best trial
+indexes = np.argsort(values)
+for i in [0]:  #choose the best-model here, e.g. [0], or [1]
+    trial = study.trials[indexes[i]]
+    print("\nTrial number {}".format(trial.number))
+    print("Value: %.5e"%trial.value)
+    print(" Params: ")
+    for key, value in trial.params.items():
+        print("    {}: {}".format(key, value))
+    n_layers = trial.params['n_layers']
+    lr       = trial.params['lr']
+    wd       = trial.params['wd']
+    hidden   = np.zeros(n_layers, dtype=np.int32)
+    dr       = np.zeros(n_layers, dtype=np.float32)
+    for i in range(n_layers):
+        hidden[i] = trial.params['n_units_l%d'%i]
+        dr[i]     = trial.params['dropout_l%d'%i]
+    fmodel = 'models/model_%d.pt'%trial.number
+
+# generate the architecture
+model = architecture.dynamic_model2(input_size, output_size, n_layers, hidden, dr)
+model.to(device)    
+
+# load best-model, if it exists
+print('Loading model...')
+if os.path.exists(fmodel):  
+    model.load_state_dict(torch.load(fmodel, map_location=torch.device(device)))
+else:  
+    raise Exception('model doesnt exists!!!')
+
+# get the data
+test_loader = data.create_dataset(mode, seed, f_Pk, f_Pk_norm, f_params, 
+                                  batch_size, shuffle=False, workers=workers)
+test_points = 0
+for x,y in test_loader:  test_points += x.shape[0]
+
+# define the matrix containing the true and predicted value of the parameters + errors
+params  = len(g)
+results = np.zeros((test_points, 3*params), dtype=np.float32)
+
+# test the model
+test_loss1, test_loss = torch.zeros(len(g)).to(device), 0.0
+test_loss2, points    = torch.zeros(len(g)).to(device), 0
+model.eval()
+for x, y in test_loader:
+    with torch.no_grad():
+        bs    = x.shape[0]         #batch size
+        x     = x.to(device)       #maps
+        y     = y.to(device)[:,g]  #parameters
+        p     = model(x)           #NN output
+        y_NN  = p[:,g]             #posterior mean
+        e_NN  = p[:,h]             #posterior std
+        loss1 = torch.mean((y_NN - y)**2,                axis=0)
+        loss2 = torch.mean(((y_NN - y)**2 - e_NN**2)**2, axis=0)
+        loss  = torch.mean(torch.log(loss1) + torch.log(loss2))
+        test_loss1 += loss1*bs
+        test_loss2 += loss2*bs
+        results[points:points+bs,0*params:1*params] = y.cpu().numpy()
+        results[points:points+bs,1*params:2*params] = y_NN.cpu().numpy()
+        results[points:points+bs,2*params:3*params] = e_NN.cpu().numpy()
+        points     += bs
+test_loss = torch.log(test_loss1/points) + torch.log(test_loss2/points)
+test_loss = torch.mean(test_loss).item()
+print('Test loss:', test_loss)
+
+# denormalize results here
+minimum = np.array([0.1, 0.6, 0.5, -1.0, -1.0, 80.])
+maximum = np.array([0.5, 1.0, 2.0,  1.0,  1.0, 1000.])
+results[:,0*params:1*params] = results[:,0*params:1*params]*(maximum-minimum)+minimum
+results[:,1*params:2*params] = results[:,1*params:2*params]*(maximum-minimum)+minimum
+results[:,2*params:3*params] = results[:,2*params:3*params]*(maximum-minimum)
+
+# save results to file
+np.savetxt(fout, results)