updates 12/10

data/signals.py

@@ -174,6 +174,18 @@ def fetch_nstx_data(signal_path, shot_num, c):
     mapping_range=(0, 1), num_channels=profile_num_channels,
     data_avail_tolerances=[0.05, 0.02])
 
+etemp_profilet = ProfileSignal(
+    "Electron temperature profile tol",
+    ["ppf/hrts/te", "ZIPFIT01/PROFILES.ETEMPFIT"], [jet, d3d],
+    mapping_paths=["ppf/hrts/rho", None], causal_shifts=[0, 10],
+    mapping_range=(0, 1), num_channels=profile_num_channels,
+    data_avail_tolerances=[0.05, 0.029])
+edens_profilet = ProfileSignal(
+    "Electron density profile tol",
+    ["ppf/hrts/ne", "ZIPFIT01/PROFILES.EDENSFIT"], [jet, d3d],
+    mapping_paths=["ppf/hrts/rho", None], causal_shifts=[0, 10],
+    mapping_range=(0, 1), num_channels=profile_num_channels,
+    data_avail_tolerances=[0.05, 0.029])
 # d3d only:
 # etemp_profile = ProfileSignal(
 #     "Electron temperature profile", ["ZIPFIT01/PROFILES.ETEMPFIT"], [d3d],
@@ -262,26 +274,45 @@ def fetch_nstx_data(signal_path, shot_num, c):
     "q95 safety factor", ['ppf/efit/q95', "EFIT01/RESULTS.AEQDSK.Q95"],
     [jet, d3d], causal_shifts=[15, 10], normalize=False,
     data_avail_tolerances=[0.03, 0.02])
+q95t = Signal(
+    "q95 safety factor tol", ['ppf/efit/q95', "EFIT01/RESULTS.AEQDSK.Q95"],
+    [jet, d3d], causal_shifts=[15, 10], normalize=False,
+    data_avail_tolerances=[0.03, 0.029])
 
 # "d3d/ipsip" was used before, ipspr15V seems to be available for a
 # superset of shots.
 ip = Signal("plasma current", ["jpf/da/c2-ipla", "d3d/ipspr15V"],
             [jet, d3d], is_ip=True)
+
+ipt = Signal("plasma current tol", ["jpf/da/c2-ipla", "d3d/ipspr15V"],
+            [jet, d3d], is_ip=True,data_avail_tolerances=[0.029, 0.029])
 iptarget = Signal("plasma current target", ["d3d/ipsiptargt"], [d3d])
+iptargett = Signal("plasma current target tol", ["d3d/ipsiptargt"], [d3d],data_avail_tolerances=[0.029])
 iperr = Signal("plasma current error", ["d3d/ipeecoil"], [d3d])
+iperrt = Signal("plasma current error tol", ["d3d/ipeecoil"], [d3d],data_avail_tolerances=[0.029])
 
 li = Signal("internal inductance", ["jpf/gs/bl-li<s", "d3d/efsli"], [jet, d3d])
+lit = Signal("internal inductance tol", ["jpf/gs/bl-li<s", "d3d/efsli"], [jet, d3d],data_avail_tolerances=[0.029, 0.029])
 lm = Signal("Locked mode amplitude", ['jpf/da/c2-loca', 'd3d/dusbradial'],
             [jet, d3d])
+lmt = Signal("Locked mode amplitude tol", ['jpf/da/c2-loca', 'd3d/dusbradial'],
+            [jet, d3d],data_avail_tolerances=[0.029, 0.029])
 dens = Signal("Plasma density", ['jpf/df/g1r-lid:003', 'd3d/dssdenest'],
               [jet, d3d], is_strictly_positive=True)
+denst = Signal("Plasma density tol", ['jpf/df/g1r-lid:003', 'd3d/dssdenest'],
+              [jet, d3d], is_strictly_positive=True,data_avail_tolerances=[0.029, 0.029])
 energy = Signal("stored energy", ['jpf/gs/bl-wmhd<s', 'd3d/efswmhd'],
                 [jet, d3d])
+energyt = Signal("stored energy tol", ['jpf/gs/bl-wmhd<s', 'd3d/efswmhd'],
+                [jet, d3d],data_avail_tolerances=[0.029, 0.029])
 # Total beam input power
 pin = Signal("Input Power (beam for d3d)", ['jpf/gs/bl-ptot<s', 'd3d/bmspinj'],
              [jet, d3d])
+pint = Signal("Input Power (beam for d3d) tol", ['jpf/gs/bl-ptot<s', 'd3d/bmspinj'],
+             [jet, d3d],data_avail_tolerances=[0.029, 0.029])
 
 pradtot = Signal("Radiated Power", ['jpf/db/b5r-ptot>out'], [jet])
+pradtott = Signal("Radiated Power tol", ['jpf/db/b5r-ptot>out'], [jet],data_avail_tolerances=[0.029])
 # pradtot = Signal("Radiated Power", ['jpf/db/b5r-ptot>out',
 # 'd3d/'+r'\prad_tot'], [jet,d3d])
 # pradcore = ChannelSignal("Radiated Power Core", [ 'd3d/' + r'\bol_l15_p']
@@ -294,17 +325,27 @@ def fetch_nstx_data(signal_path, shot_num, c):
 pradedge = ChannelSignal("Radiated Power Edge",
                          ['ppf/bolo/kb5h/channel10', 'd3d/' + r'\bol_l03_p'],
                          [jet, d3d])
+pradcoret = ChannelSignal("Radiated Power Core tol",
+                         ['ppf/bolo/kb5h/channel14', 'd3d/' + r'\bol_l15_p'],
+                         [jet, d3d],data_avail_tolerances=[0.029, 0.029])
+pradedget = ChannelSignal("Radiated Power Edge tol" ,
+                         ['ppf/bolo/kb5h/channel10', 'd3d/' + r'\bol_l03_p'],
+                         [jet, d3d],data_avail_tolerances=[0.029, 0.029])
 # pechin = Signal("ECH input power, not always on", ['d3d/pcechpwrf'], [d3d])
 pechin = Signal("ECH input power, not always on",
                 ['RF/ECH.TOTAL.ECHPWRC'], [d3d])
+pechint = Signal("ECH input power, not always on tol",
+                ['RF/ECH.TOTAL.ECHPWRC'], [d3d],data_avail_tolerances=[0.029])
 
 # betan = Signal("Normalized Beta", ['jpf/gs/bl-bndia<s', 'd3d/efsbetan'],
 # [jet, d3d])
 betan = Signal("Normalized Beta", ['d3d/efsbetan'], [d3d])
+betant = Signal("Normalized Beta tol", ['d3d/efsbetan'], [d3d],data_avail_tolerances=[0.029])
 energydt = Signal(
     "stored energy time derivative", ['jpf/gs/bl-fdwdt<s'], [jet])
 
 torquein = Signal("Input Beam Torque", ['d3d/bmstinj'], [d3d])
+torqueint = Signal("Input Beam Torque tol", ['d3d/bmstinj'], [d3d],data_avail_tolerances=[0.029])
 tmamp1 = Signal("Tearing Mode amplitude (rotating 2/1)", ['d3d/nssampn1l'],
                 [d3d])
 tmamp2 = Signal("Tearing Mode amplitude (rotating 3/2)", ['d3d/nssampn2l'],
@@ -314,27 +355,26 @@ def fetch_nstx_data(signal_path, shot_num, c):
 tmfreq2 = Signal("Tearing Mode frequency (rotating 3/2)", ['d3d/nssfrqn2l'],
                  [d3d])
 ipdirect = Signal("plasma current direction", ["d3d/iptdirect"], [d3d])
+ipdirectt = Signal("plasma current direction tol", ["d3d/iptdirect"], [d3d],data_avail_tolerances=[0.029])
 
 # for downloading, modify this to preprocess shots with only a subset of
 # signals. This may produce more shots
 # since only those shots that contain all_signals contained here are used.
-# all_signals = {'q95':q95, 'li':li, 'ip':ip,
-# 'betan':betan, 'energy':energy, 'lm':lm, 'dens':dens,
-# 'pradcore':pradcore, 'pradedge':pradedge, 'pradtot':pradtot,
-# 'pin':pin, 'torquein':torquein,
-# 'tmamp1':tmamp1, 'tmamp2':tmamp2, 'tmfreq1':tmfreq1, 'tmfreq2':tmfreq2,
-# 'pechin':pechin, 'energydt':energydt,
-# 'etemp_profile':etemp_profile, 'edens_profile':edens_profile,
-# 'ipdirect':ipdirect, 'iptarget':iptarget, 'iperr':iperr,
-# 'etemp_profile_spatial':etemp_profile_spatial,
-# 'edens_profile_spatial':edens_profile_spatial,
-
-# 'etemp':etemp
-# }
 
 # Restricted subset to those signals that are present for most shots. The
 # idea is to remove signals that cause many shots to be dropped from the
 # dataset.
+all_signals_gar18 = {
+    'q95t': q95t, 'lit': lit, 'ipt': ipt, 'betant': betant, 'energyt': energyt, 'lmt': lmt,
+    'denst': denst, 'pradcoret': pradcoret,
+    'pradedget': pradedget, 'pint': pint,
+    'torqueint': torqueint,
+    'ipdirectt': ipdirectt, 'iptargett': iptargett,
+    'iperrt': iperrt,
+    'etemp_profilet': etemp_profilet, 'edens_profilet': edens_profilet,
+}
+
+
 all_signals = {
     'q95': q95, 'li': li, 'ip': ip, 'betan': betan, 'energy': energy, 'lm': lm,
     'dens': dens, 'pradcore': pradcore,
@@ -353,6 +393,9 @@ def fetch_nstx_data(signal_path, shot_num, c):
     # 'q_psi_profile':q_psi_profile}
 }
 
+
+
+
 # new signals are not downloaded yet
 
 
@@ -362,9 +405,10 @@ def fetch_nstx_data(signal_path, shot_num, c):
 
 all_signals_restricted = all_signals
 
-g.print_unique('All signals (determines which signals are downloaded'
-               ' & preprocessed):')
-g.print_unique(all_signals.values())
+#g.print_unique('All signals (determines which signals are downloaded'
+#               ' & preprocessed):')
+#g.print_unique(all_signals.values())
+
 
 fully_defined_signals = {
     sig_name: sig for (sig_name, sig) in all_signals_restricted.items() if (
@@ -382,6 +426,10 @@ def fetch_nstx_data(signal_path, shot_num, c):
     sig_name: sig for (sig_name, sig) in all_signals_restricted.items() if (
         sig.is_defined_on_machine(d3d))
 }
+d3d_signals_gar18 = {
+    sig_name: sig for (sig_name, sig) in all_signals_gar18.items() if (
+        sig.is_defined_on_machine(d3d))
+}
 d3d_signals_0D = {
     sig_name: sig for (sig_name, sig) in all_signals_restricted.items() if (
         (sig.is_defined_on_machine(d3d) and sig.num_channels == 1))

examples/conf.yaml

@@ -10,7 +10,7 @@ paths:
     signal_prepath: '/signal_data/' # /signal_data/jet/
     shot_list_dir: '/shot_lists/'
     tensorboard_save_path: '/Graph/'
-    data: d3d_data_0D # 'd3d_to_jet_data' # 'd3d_to_jet_data' #  'jet_to_d3d_data' # jet_data
+    data: d3d_data ##_0D # 'd3d_to_jet_data' # 'd3d_to_jet_data' #  'jet_to_d3d_data' # jet_data
     # if specific_signals: [] left empty, it will use all valid signals defined on a machine. Only use if need a custom set
     specific_signals: [] # ['q95','li','ip','betan','energy','lm','pradcore','pradedge','pradtot','pin','torquein','tmamp1','tmamp2','tmfreq1','tmfreq2','pechin','energydt','ipdirect','etemp_profile','edens_profile']
     executable: "mpi_learn.py"
@@ -74,23 +74,23 @@ model:
         mlp_regularization: 0.0001
         skip_train: False # should a finished model be loaded if available
     # length of LSTM memory
-    pred_length: 200
+    pred_length: 256
     pred_batch_size: 128
     # TODO(KGF): optimize length of LSTM memory
-    length: 128
+    length: 256
     skip: 1
     # hidden layer size
     # TODO(KGF): optimize size of RNN layers
     rnn_size: 200
     # size 100 slight overfitting, size 20 no overfitting. 200 is not better than 100. Prediction much better with size 100, size 20 cannot capture the data.
     rnn_type: 'LSTM'
     # TODO(KGF): optimize number of RNN layers
-    rnn_layers: 2
+    rnn_layers: 1
     num_conv_filters: 128
     size_conv_filters: 3
     num_conv_layers: 3
     pool_size: 2
-    dense_size: 128
+    dense_size: 64
     extra_dense_input: False
     # have not found a difference yet
     optimizer: 'adam'
@@ -103,7 +103,7 @@ model:
     lr: 0.00002 # 0.00001 # 0.0005 # for adam plots 0.0000001 # 0.00005 # 0.00005 # 0.00005
     lr_decay: 0.97 # 0.98 # 0.9
     stateful: True
-    return_sequences: True
+    return_sequences: True # True
     dropout_prob: 0.1
     # only relevant if we want to do MPI training. The number of steps with a single replica
     warmup_steps: 0

examples/slurm.cmd

@@ -11,7 +11,7 @@
 # Each node = 2.4 GHz Xeon Broadwell E5-2680 v4 + 4x 1328 MHz P100 GPU
 
 module load anaconda3
-conda activate my_env
+conda activate pplori
 module load cudatoolkit
 module load cudnn
 module load openmpi/cuda-8.0/intel-17.0/3.0.0/64

plasma/conf_parser.py

@@ -27,12 +27,21 @@ def parameters(input_file):
         base_path = output_path
 
         params['paths']['base_path'] = base_path
-        params['paths']['signal_prepath'] = (
-            base_path + params['paths']['signal_prepath'])
+        if isinstance(params['paths']['signal_prepath'],list):
+            print('reading from multiple data folder!**********************************************')
+
+            params['paths']['signal_prepath'] = [base_path+s for s in params['paths']['signal_prepath']]
+        else:
+            params['paths']['signal_prepath']=base_path+params['paths']['signal_prepath']
         params['paths']['shot_list_dir'] = (
             base_path + params['paths']['shot_list_dir'])
         params['paths']['output_path'] = output_path
-        h = myhash_signals(sig.all_signals.values())
+        if params['paths']['data']=='d3d_data_gar18':
+           h = myhash_signals(sig.all_signals_gar18.values())
+        elif params['paths']['data']=='d3d_data_garbage':
+           h = myhash_signals(sig.all_signals_gar18.values())*2
+        else:   
+           h = myhash_signals(sig.all_signals.values())
         params['paths']['global_normalizer_path'] = (
             output_path
             + '/normalization/normalization_signal_group_{}.npz'.format(h))
@@ -85,7 +94,11 @@ def parameters(input_file):
         # params['model']['loss'] = params['data']['target'].loss
 
         # signals
-        params['paths']['all_signals_dict'] = sig.all_signals
+        if params['paths']['data'] in ['d3d_data_gar18','d3d_data_garbage']:
+           params['paths']['all_signals_dict'] = sig.all_signals_gar18
+        else:
+           params['paths']['all_signals_dict'] = sig.all_signals
+
         # assert order
         # q95, li, ip, lm, betan, energy, dens, pradcore, pradedge, pin,
         # pechin, torquein, ipdirect, etemp_profile, edens_profile
@@ -136,6 +149,7 @@ def parameters(input_file):
             sig.d3d, params['paths']['shot_list_dir'],
             ['d3d_clear_data_avail.txt', 'd3d_disrupt_data_avail.txt'],
             'd3d data since shot 125500')
+        d3d_full_new = ShotListFiles(sig.d3d,params['paths']['shot_list_dir'],['shots_since_2016_clear.txt','shots_since_2016_disrupt.txt'],'d3d data since shot 125500')
         d3d_jenkins = ShotListFiles(
             sig.d3d, params['paths']['shot_list_dir'],
             ['jenkins_d3d_clear.txt', 'jenkins_d3d_disrupt.txt'],
@@ -224,6 +238,49 @@ def parameters(input_file):
                 'etemp_profile': sig.etemp_profile,
                 'edens_profile': sig.edens_profile,
             }
+        elif params['paths']['data'] in ['d3d_data_gar18','d3d_data_garbage']:
+            params['paths']['shot_files'] = [d3d_full_new]
+            params['paths']['shot_files_test'] = []
+            params['paths']['use_signals_dict'] = {
+                'q95t': sig.q95t,
+                'lit': sig.lit,
+                'ipt': sig.ipt,
+                'lmt': sig.lmt,
+                'betant': sig.betant,
+                'energyt': sig.energyt,
+                'denst': sig.denst,
+                'pradcoret': sig.pradcoret,
+                'pradedget': sig.pradedget,
+                'pint': sig.pint,
+                'torqueint': sig.torqueint,
+                'ipdirectt': sig.ipdirectt,
+                'iptargett': sig.iptargett,
+                'iperrt': sig.iperrt,
+                'etemp_profilet': sig.etemp_profilet,
+                'edens_profilet': sig.edens_profilet,
+            }
+
+        elif params['paths']['data'] == 'd3d_data_new':
+            params['paths']['shot_files'] = [d3d_full_new]
+            params['paths']['shot_files_test'] = [] 
+            params['paths']['use_signals_dict'] = {
+                'q95': sig.q95,
+                'li': sig.li,
+                'ip': sig.ip,
+                'lm': sig.lm,
+                'betan': sig.betan,
+                'energy': sig.energy,
+                'dens': sig.dens,
+                'pradcore': sig.pradcore,
+                'pradedge': sig.pradedge,
+                'pin': sig.pin,
+                'torquein': sig.torquein,
+                'ipdirect': sig.ipdirect,
+                'iptarget': sig.iptarget,
+                'iperr': sig.iperr,
+                'etemp_profile': sig.etemp_profile,
+                'edens_profile': sig.edens_profile,
+            }
         elif params['paths']['data'] == 'd3d_data_1D':
             params['paths']['shot_files'] = [d3d_full]
             params['paths']['shot_files_test'] = []

plasma/models/builder.py

@@ -25,6 +25,7 @@
 from copy import deepcopy
 from plasma.utils.downloading import makedirs_process_safe
 from plasma.utils.hashing import general_object_hash
+from plasma.models.tcn import TCN
 
 # Synchronize 2x stderr msg from TensorFlow initialization via Keras backend
 # "Succesfully opened dynamic library... libcudart" "Using TensorFlow backend."
@@ -159,8 +160,12 @@ def slicer_output_shape(input_shape, indices):
             # slicer_output_shape(s,indices_0d))(pre_rnn_input)
             pre_rnn_1D = Reshape((num_1D, len(indices_1d)//num_1D))(pre_rnn_1D)
             pre_rnn_1D = Permute((2, 1))(pre_rnn_1D)
-
-            for i in range(model_conf['num_conv_layers']):
+            if 'simple_conv' in model_conf.keys() and model_conf['simple_conv']==True:
+              for i in range(model_conf['num_conv_layers']):
+                pre_rnn_1D = Convolution1D(num_conv_filters,size_conv_filters,padding='valid',activation='relu') (pre_rnn_1D)
+                pre_rnn_1D = MaxPooling1D(pool_size) (pre_rnn_1D) 
+            else:
+              for i in range(model_conf['num_conv_layers']):
                 div_fac = 2**i
                 '''The first conv layer learns `num_conv_filters//div_fac`
                 filters (aka kernels), each of size
@@ -254,10 +259,48 @@ def slicer_output_shape(input_shape, indices):
                 activity_regularizer=l2(dense_regularization))(pre_rnn)
 
         pre_rnn_model = Model(inputs=pre_rnn_input, outputs=pre_rnn)
+        from mpi4py import MPI
+        comm = MPI.COMM_WORLD
+        task_index = comm.Get_rank()
+        if not predict and task_index==0 :
+          print('Printingout pre_rnn model.........')
+          fr=open('model_architecture.log','w')
+          ori=sys.stdout
+          sys.stdout=fr
+          pre_rnn_model.summary()
+          sys.stdout=ori
+          fr.close()        
         # pre_rnn_model.summary()
         x_input = Input(batch_shape=batch_input_shape)
-        x_in = TimeDistributed(pre_rnn_model)(x_input)
-        for _ in range(model_conf['rnn_layers']):
+        if num_1D>0 or (
+                'extra_dense_input' in model_conf.keys()
+                and model_conf['extra_dense_input']):
+            x_in = TimeDistributed(pre_rnn_model)(x_input)
+        else:
+            x_in=x_input
+
+###################TCN model##################################
+        if 'keras_tcn' in model_conf.keys() and model_conf['keras_tcn']==True:
+          print('Building TCN model....')
+          tcn_layers=model_conf['tcn_layers']
+          tcn_dropout=model_conf['tcn_dropout']
+          nb_filters=model_conf['tcn_hidden']
+          kernel_size=model_conf['kernel_size_temporal']
+          nb_stacks=model_conf['tcn_nbstacks']
+          use_skip_connections=model_conf['tcn_skip_connect']
+          activation=model_conf['tcn_activation']
+          use_batch_norm=model_conf['tcn_batch_norm']
+          for _ in range(model_conf['tcn_pack_layers']):
+            x_in=TCN(
+               use_batch_norm=use_batch_norm,activation=activation,
+               use_skip_connections=use_skip_connections,
+               nb_stacks=nb_stacks,kernel_size=kernel_size,
+               nb_filters=nb_filters,num_layers=tcn_layers,
+               dropout_rate=tcn_dropout)(x_in)
+            x_in = Dropout(dropout_prob) (x_in)
+        else: 
+###################TCN model##################################
+          for _ in range(model_conf['rnn_layers']):
             x_in = rnn_model(
                 rnn_size, return_sequences=return_sequences,
                 # batch_input_shape=batch_input_shape,