#!/bin/bash
# If sourced, will only load its functions, otherwise will run the setup

################################################################################
# SPECIFC slurm
#SBATCH --job-name=htexplo
#SBATCH --ntasks=32               # nb process mpi
#SBATCH --ntasks-per-node=32      # mpi / node
#SBATCH --cpus-per-task=1         # omp/mpi
#SBATCH --hint=nomultithread      # 1 thread /core (no hyperthreading)
#SBATCH --output=outhtexplo%j     # output file
#SBATCH --error=outhtexplo%j      # error file
#SBATCH --time=04:00:00
################################################################################

ulimit -s unlimited
set -eu  # exit on most (not all!) error codes, and error on unset variables
export LC_ALL=C # set locale for computation with awk

################################################################################
# HighTune bench.sh
#
# Author : Hourdin et al., frederic.hourdin@lmd.ipsl.fr
# Started in 2017
# Originally written as a bench of the htexplo tools, it has evolved
# in an an automatic tool that can run most of the applications.
# It can serve also as a guide on how to use it.
#
################################################################################
#
# General structure and main functions called
#
# 0) Initialisations
#
# for N in 1 ... Nwave ; do
#
#       2) if (N==1) param2R.sh -> ModelParam.R
#          if (N==1)  -> Waves[N].RData +  Par1D_Wave[N].asc
#          else          Waves[N].RData -> Par1D_Wave[N].asc
#
#       3) serie_[MODEL].sh
#             Par1D_Wave3.asc -> Simulations
#             Results on WAVE[N]/[CASE]/[SUBCASE]/SCM*nc
#
#       4) compute_metrics.sh
#             simulations -> metrics csv
#          htune_csv2Rdata.R
#             csv to Rdata including adding old simulations.
#
#       5) Emulator + history matching
#            htune_Emulating_Multi_Metric_Multi_LHS.R
#            Emulators stored in :
#            -> WAVE[N]/EMULATOR_MULT_METRIC_wave[N].RData
#            -> WAVE[N]/EMULATOR_MULT_METRIC_wave[N]_mogp
#            New sample of parameter vectors in NROY[N]
#            -> Waves[N+1].RData
#
# done
#
################################################################################
#  The htexplo tools are inherited from the HighTune project,
#  and are described in a 3-part paper titled
#  "Process-based climate model development harnessing machine learning"
#  Part I : General presentation
#  https://www.lmd.jussieu.fr/~hourdin/PUBLIS/ItuneI2020.pdf
#  Part II : An application to the LMDZ model
#  https://www.lmd.jussieu.fr/~hourdin/PUBLIS/Hourdin_HighTune_2021_Proof.pdf
#  Part III : An application to raditive transfer comupations
#
#  Results of part II can be reproduced by running the following commands :
#  Section 4 :
# ./bench.sh -wdir ArtII1 -waves "`seq 1 20`" -param param_ArtII1 -metrics "ARMCU_REF_zav-400-600-theta_7_9,ARMCU_REF_zav-400-600-qv_7_9,ARMCU_REF_nebmax_7_9,RICO_REF_nebmax_19_25,SANDU_REF_neb4zave_50_60"
#  Section 5 :
#  ./bench.sh -wdir ArtII2 -waves "`seq 1 40`" -param param_ArtII2 -sample_size 3000000 -sample_size_next_design 90 -metrics "ARMCU_REF_zav-400-600-theta_7_9,ARMCU_REF_zav-400-600-qv_7_9,IHOP_REF_zav-400-600-theta_7_9,ARMCU_REF_nebzave_7_9,ARMCU_REF_neb4zave_7_9,ARMCU_REF_nebmax_7_9,RICO_REF_nebmax_19_25,SANDU_SLOW_neb4zave_50_60,SANDU_REF_neb4zave_50_60,SANDU_REF_nebzave_50_60,SANDU_FAST_neb4zave_50_60"

################################################################################

function read_arguments_and_defaults() {
  echo '====================================================================='
  echo "O/ default values and reading arguments"
  echo '====================================================================='


  # 0.1/ Default values
  wdir=""
  expdir="WORK"
  installdir=$(readlink -f "../install_shared")

  model="LMDZ"
  metrics="RICO_REF_nebmax_9_9"
  param="param"
  waves="1" # could be waves=$(seq 1 15), waves="1 2 3"
  serie=""
  sample_size=300000
  sample_size_next_design=45
  sepm="_"
  cutoff_val=(3 2.5 2)
  cutoff_wav=(1  5  8)

  GCM="no"
  # Managing GCM GCM waves
  # GCM = no : SCM/LES only
  # GCM = pre : preparation of a wave. From design to metrics computation in 1D
  #            before breaking for running GCM, computing associated metrics and combining them with 1D
  # GCM = post : running emulator and preparing the next wave (pre)

  # to only print commands and not execute them
  dryrun=0
  do_plot=1

  opt_model=""
  metrics_file=""
  seed=""
  tau=0

  # 0.3/ options
  while (($# > 0)) ; do
    case $1 in
      -installdir) installdir=$2; shift 2 ;;
      -expdir) expdir=$2; shift 2 ;;
      -tau) tau=$2; shift 2 ;;
      -serie) serie=$2 ; shift 2 ;; # useful for ecRad runs
      -wdir) wdir=$2 ; shift 2 ;;
      -param) param=$2 ; shift 2 ;;
      -sample_size|-ss) sample_size=$2 ; shift 2 ;;
      -sample_size_next_design|-ssnd) sample_size_next_design=$2 ; shift 2 ;;
      -waves)  waves="$2"  ; shift 2 ;;
      -GCM)  GCM="$2"  ; shift 2 ;;
      -model) model=$2 ; shift 2 ;;
      -cutoff_val|-cv) IFS=" " read -r -a cutoff_val <<< "$2"; shift 2 ;;
      -cutoff_wav|-cw) IFS=" " read -r -a cutoff_wav <<< "$2" ; shift 2 ;;
      -opt_model) opt_model="$2" ; shift 2 ;;
      -metrics|-m) metrics="${2//,/ }" ; shift 2 ;;
      -metrics_file|-mf) metrics_file=$2 ; shift 2
        if [[ ! -f $metrics_file ]] ; then echo "$metrics_file does not exists" ; exit 1 ; fi
        metrics="$( grep ${sepm}'.*'${sepm}'.*'${sepm}'.*'${sepm}'.*,.*,' "$metrics_file" | cut -d, -f1 )"
        echo "Metrics read in $metrics_file : $metrics" ;;
      -seed) seed="-seed $2" ; shift 2 ;;
      -dry) dryrun=1 ; shift ;;
      -noplot) do_plot=0  ; shift ;;
      -h|-help|--help) echo "Usage: $0 [-param param_file] [-waves \"1 [2 3 ...]\"] [-wdir DIRNAME] [-sample_size sample_size] [-model model] [-metrics metrics1,metrics2,...] [-noplot] or directly $0 model" ; cat <<eod
-param param_file : param_file contains the name, the min/max/nominal values, and the mode of exploration Linear/Log
            of the parameters (default <$param>)
-wdir WDIR        : the history matching sequence will be run on \$EXPDIR/\$WDIR (default <$wdir>)
-waves WAVES      : WAVES is a sequence of numbers. 1 ; "1 2 3" ; "\`seq 1 20\`" (default <$waves>)
            Can start at N+1 if waves 1 to N are already done
-sample_size|-ss SAMPLESIZE : sample size for the NROY graphics (default <$sample_size>)
-sample_size_next_design|-ssnd SAMPLESIZENEX : sample size for next design (default <$sample_size_next_design>)
-model MODEL      : name of MODEL, available on models/  (default <$model>)
-cutoff_val|-cv   : list of cutoff values (default <${cutoff_val[*]}>)
-cutoff_wav|-cw   : list of waves associated with cutoff value changes  (default <${cutoff_wav[*]}>)
-opt_model opt    : option to be given to serie_$model (default <$opt_model>)
-metrics|-m METRICS  : METRICS is a list of metrics separated by "," or " " (default <$metrics>)
-metrics_file|-mf METRICS  : METRICS is a list of metrics separated by "," or " " (default <$metrics_file>)
-seed N           : imposing the seed for the random number to ensure reproductibility. Use it for quality control only !  (default <$seed>)
-installdir       : Where to install the tuning environment (default <$installdir>)
eod
        exit 0 ;;
      *) model=$1 ; shift ;;
    esac
  done

  if [[ $wdir = "" ]] ; then wdir="BENCH$model" ; fi

  #list of available 1D cases
  list_case="ARMCU/REF ARMCU/MESONH ARMCU/MNHRAD ARMCU/MNHRADSURF RICO/REF RICO/MESONH FIRE/REF BOMEX/REF AYOTTE/24SC AYOTTE/05WC IHOP/REF SANDU/SLOW SANDU/REF AMMA/REF SANDU/FAST RCEOCE/REF GABLS1/REF GABLS4/STAGE3 GABLS4/STAGE3SHORT CINDY/REF EUROCS/REF KB2006/MESONH KB2006/REF MPACE/REF ISDAC/REF COMBLE/REF"

  if [[ ! -f models/$model/$param ]] ; then echo "No $param file found" ; exit 1 ; fi

  # controling cutoff consistency
  ncutoff=${#cutoff_val[@]}
  if [[ $ncutoff != "${#cutoff_wav[@]}" ]]; then echo "Pb in cutoff specification"
    echo "vals ${cutoff_val[*]}" ; echo "waves ${cutoff_wav[*]}" ; exit 1
  fi
  cutoff_wav[0]=1

}

function run() {
  # the run function will print command and exe or not depending on value of $dryrun
  echo "$@"
  if [[ $dryrun == 0 ]] ; then
    "$@"
  fi
}

function run_post() {
  echo '====================================================================='
  echo "1A/ Continuing an experiment after 3D GCM. Nothing to initialize"
  echo '====================================================================='
  echo "CAS POST GCM $waves"
  if [[ $(echo "$waves" | wc -w) != "1" ]] ; then
    echo "Just one wave when GCM is post" ; exit 1
  fi
}

function run_resume() {
  echo '====================================================================='
  echo "1B/ Continuing an existing experiment"
  echo '====================================================================='

  source "$installdir/env.sh"
  echo "Directory WORK/$wdir already exists. Do you want to continue (Y/n) ?"
  local answ
  read -r answ
  #echo "Forcing answer to Y"; answ="Y"  # TODO (Amaury) revert to read and when called use a piped "Y" answer
  if [[ $answ = "n" ]] ; then
      exit 0
  fi

  cd "$expdir/$wdir"
  local first_wtbd=$(echo $waves | awk ' { print $1 } ')
  # Checking whether the sampling is available for WAVE$first_wtbd
  if [[ ! -f WAVE$first_wtbd/Wave${first_wtbd}.RData_orig \
    && ! -f Wave${first_wtbd}.RData ]] ; then
    echo "No Wave${first_wtbd}.RData file available" ; exit 1
  fi

  local wavesdone lastwdone
  wavesdone="" ;
  for w in $(seq 1 100) ; do
    if [[ -d WAVE$w ]] ; then wavesdone="$wavesdone WAVE$w" ; fi
  done
  # Case were some of the waves tbd are there already
  lastwdone=$(echo "$wavesdone" | awk ' { print $NF } ' | sed -e 's/WAVE//')
  echo first_wtbd $first_wtbd
  echo wavesdone $wavesdone
  echo lastwdone $lastwdone
  if [[ "$lastwdone" -ne "$(( first_wtbd - 1 ))" ]] ; then
    echo "WAVES $wavesdone were already present"
    echo "From which wave do you really like to restart, the following\
               ones being be saved as WAVEN_$$ (0 to stop)"
    read -r first_wtbd
    # waves="$(seq "$first_wtbd" "$(echo "$waves" | awk ' { print $NF } ')")"
    waves=$(seq $first_wtbd $(echo $waves | awk ' { print $NF } '))

  fi

  if [[ ! -f Wave${first_wtbd}.RData ]] ; then
    mv "WAVE${first_wtbd}/Wave${first_wtbd}.RData_orig" "Wave${first_wtbd}.RData"
  fi

  echo '--------------------------------------------------------------------'
  echo "1A.1 saving previously run waves"
  for w in $waves ; do
    if [[ -d WAVE$w ]] ; then
      echo "saving wave $w in wave ${w}_$$"
      mv "WAVE$w" "WAVE${w}_$$"
    fi
  done # Saving previously run waves
  cd -
}

#-------------------------------------------------------------------------------
function setup_new_experiment() {
#-------------------------------------------------------------------------------
  echo '======================================================================'
  echo "1C/ Running setup.sh $model $wdir for a new experiment"
  echo '======================================================================'

  # 0.4/ Temporary trick to be able to compare old (ExeterUQ) and new
  #     (ExeterUQ_MOGP) Exeter tools
  #     To be removed as soon as the new versions are stabilized (2020/07/20)
  ExeterUQ="ExeterUQ_MOGP" # ExeterUQ=ExeterUQ to come back to the old version
  sed -e 's/^ExeterUQ=.*.$/ExeterUQ='$ExeterUQ'/' "setup.sh" > setup.sh.tmpsed
  mv setup.sh.tmpsed setup.sh
  chmod +x setup.sh

  mkdir -p log

  local models_to_install="$model"
  set +e
  echo $metrics|sed -e "s/ /\n/g"|awk -F"_" '{print $1}'|grep -w RAD > /dev/null 2>&1
  is_metric_rad=$(( 1-$? ))
  set -e
  echo $is_metric_rad
  if [[ $is_metric_rad -eq 1 && $model != "ECRAD" ]] ; then
    models_to_install="$models_to_install ECRAD"
  fi
  echo $models_to_install

  local model_to_install
  for model_to_install in $models_to_install; do
    local setup_cmd
    setup_cmd="./setup.sh -installdir $installdir -model $model_to_install -exp $wdir -expdir $expdir &> $(pwd)/log/setup_$model_to_install.log$$"
    echo "Running $setup_cmd"
    if ! eval "$setup_cmd" ; then
      echo "Error in setup, tail log/setup_$model_to_install.log$$ : $(tail "$(pwd)/log/setup_$model_to_install.log$$")" ; exit 1
    fi
  done

  source "$installdir/env.sh"

  if [[ "$metrics_file" != "" ]] ; then cp "$metrics_file" "$expdir/$wdir/metrics_def.csv" ; fi
}

#-------------------------------------------------------------------------------
function build_wave_design() {
#-------------------------------------------------------------------------------
  echo '======================================================================'
  echo "2/ Building design for wave $wave"
  echo '======================================================================'

  if [[ $wave = 1 ]] ; then
    # Generating ModelParam.R from $param file
    run "./param2R.sh" "$param" >> $(pwd)/log/sample_$wave.log$$
    local nsample=1
    # For the first wave, the sampling is done a $nsample sub-sampling
    # of size $subsample_size with a Latin Hypercube sampling
    subsample_size=$(( $sample_size_next_design / $nsample + 1 ))
    run Rscript htune_convertDesign.R -LHCSIZE $subsample_size -NLHC "$nsample" $seed >> $(pwd)/log/sample_$wave.log$$
  else
    cp "Wave$wave.RData" "Wave$wave.RData_orig"
    echo "2B/ Computing the LHS for wave $wave"
    run Rscript htune_convertDesign.R -wave "$wave" $seed >> $(pwd)/log/sample_$wave.log$$
  fi
  mkdir "WAVE${wave}"
  mv "Wave${wave}.RDat"* "WAVE${wave}/"
  mv "Par1D_Wave${wave}.asc" "WAVE${wave}/"
}

#-------------------------------------------------------------------------------
function run_SCM_for_wave() {
#-------------------------------------------------------------------------------
  echo '======================================================================'
  echo "3/ Running the required SCM simulations for wave $wave"
  echo '======================================================================'

  # extracting the list of required simulations from the list of
  # metrics ; simus = CASE/SUBCASE

  local simus_ case_ simus_rad simus
  simus_=$( for m in $metrics ; do
    if [[ ${m:0:3} != RAD ]] ; then
      echo "$m" | awk -F"$sepm" ' {print $1"/"$2 }'
    else # RAD
      str=$(echo "$m" | awk -F"$sepm" ' {print $2"/"$3 } ')
      echo "${str:0:-3}" # remove time to keep only eg ARMCU/REF
    fi
    done | sort | uniq )

  # we only run cases available in liste_case
  simus=""
  for case_ in $simus_ ; do
    if (echo "$list_case" | grep -q "$case_") ; then
      simus="$simus $case_"
    fi
  done

  sed_i "s/WAVEN=.*/WAVEN=$wave/g" expe_setup.R

  # extracting the list of required simulations from the list of
  # metrics ; only radiative metrics (first 3 chars = RAD) :
  # simus_rad = RAD/CASE/SUBCASE
  simus_rad=$(for m in $metrics ; do if [[ ${m:0:3} == RAD ]] ; then echo "$m" | awk -F$sepm ' {print "RAD/"$2"/"$3 } ' ; fi ; done | sort | uniq)
  simus=$(echo "$simus" | xargs)  # remove extra spaces & \n
  simus_rad=$(echo "$simus_rad" | xargs)  # remove extra spaces & \n
  echo "LIST SIMU SCM: <$simus>"
  echo "LIST SIMU RAD: <$simus_rad>"

  local wave_cmd
  if [[ "$model" != ECRAD ]]; then # run scm on $simus (list of CASE/SUBCASE)
    wave_cmd="./serie_$model$serie.sh $opt_model $simus $wave &> $(pwd)/log/serie_$wave.log$$"
    if ! run eval "$wave_cmd"; then
      echo "Error during serie_$model$serie.sh"
      tail "$(pwd)/log/serie_$wave.log$$"
      exit 1
    fi
  fi

  if [[ "$simus_rad" != "" ]]; then # run ecrad on $simus_rad on 1D LES if model=ECRAD or scm runs otherwise
    wave_cmd="./serie_ECRAD$serie.sh $simus_rad $wave $model &> $(pwd)/log/serie_${wave}_ECRAD.log$$"
    if ! run eval "$wave_cmd"; then
     echo "Error during serie_ECRAD${serie}.sh"
     tail "$(pwd)/log/serie_${wave}_ECRAD.log$$"
     exit 1
    fi
  fi
}

function compute_metrics() {
  echo '===================================================================='
  echo "4/ Computing metrics running for wave $wave"
  echo '======================================================================'

  local metrics_ case_
  metrics_=""
  for m in $metrics ; do
    if [[ ${m:0:3} != RAD ]] ; then
      case_=$(echo "$m" | awk -F$sepm  ' {print $1"/"$2 }' )
    else # RAD*
      str=$(echo "$m" | awk -F$sepm  ' {print $2"/"$3 } ')
      case_=${str:0:-3}
    fi

    if echo "$list_case" | grep -q "$case_" ; then
      metrics_="$metrics_ $m"
    fi
  done

  local metrics_cmd
  metrics_cmd="./compute_metrics.sh $metrics_ -wave $wave &> $(pwd)/log/metrics_$wave.log$$"
  if ! run eval "time $metrics_cmd"; then
    echo "Error during compute_metrics.sh"
    exit 1
  fi
}

function launch_postproc() {
  local pproc_cmd
  pproc_cmd="bash post_processing_1D.sh $wave &> $(pwd)/log/out_post_processing_1D_$wave.log$$"
  run eval "$pproc_cmd" &
}

function update_sample_size() {
  echo '--------------------------------------------------------------------'
  echo "5.1/ Increasing sample_size if Nroy too small for wave > 1"
  echo '--------------------------------------------------------------------'

  function comp_sample_size() {
    if [[ $wave -gt 1  ]] ; then
      local anticipated_reduction prev_remaining_space nroy_size is_es n1 sample_size_next_min
      anticipated_reduction=1.5
      prev_remaining_space="Remaining_space_after_wave_$(( $wave - 1 )).txt"
      echo      "By $anticipated_reduction from one wave to the other"
      echo      "Reading the previous NROY fraction in $prev_remaining_space"
      if [[ ! -f $prev_remaining_space ]] ; then
        echo "Error: $prev_remaining_space does not exist" ; exit 1
      else
        nroy_size=$(cat $prev_remaining_space)
        if [[ "$nroy_size" = 0 || "$nroy_size" = "" ]] ; then
          echo "Error: $prev_remaining_space empty or containing 0" ; exit 1
        else
          # bc - l doesn't handle scientific notation :
          is_es=$(echo "$nroy_size" | grep e || true)
          if [[ -n $is_es ]] ; then
            #convert nroy_size to be handled by bc -l
            n1=$(echo "$nroy_size" | sed 's/e/\\*10\\^/' | sed 's/+//')
          else
            n1=$nroy_size
          fi

          nroy_size_next=$(echo "$n1 / $anticipated_reduction" | bc -l)
          sample_size_next_min=$(echo "$sample_size_next_design / $nroy_size_next" | bc -l | cut -d. -f1)
          if [[ $sample_size -lt $sample_size_next_min ]] ; then
            sample_size=$sample_size_next_min
          fi
          echo "New sample_size: $sample_size"
        fi
      fi
    fi
  }
  run comp_sample_size
}

function use_previous_waves_results() {
  echo '----------------------------------------------------------------------'
  echo "5.2 Introducing results of the previous best simulations"
  echo '----------------------------------------------------------------------'

  cmd='cat "WAVE'$wave'/Par"*asc > "WAVE'$wave'/Params.asc"'
  echo $cmd
  eval $cmd

  cmd='cat "WAVE'$wave'/metrics_WAVE"*csv > "WAVE'$wave'/Metrics.csv"'
  echo $cmd
  eval $cmd

  local include_previous_waves="None"
  local iw sim

  case $include_previous_waves in
    None)  echo "Using only simulations from last wave" ;;
    All)   for iw in $(seq 1 $(( $wave - 1 ))) ; do
        tail -n +2 "WAVE${iw}/Par1D_Wave${iw}.asc" >> "WAVE${wave}/Params.asc"
        tail -n +2 "WAVE${iw}/metrics_WAVE${iw}_${iw}.csv" >> "WAVE${wave}/Metrics.csv"
      done ;;
    Bests) if [[ $wave -gt 1 && -f PrevSimulationsToBeIncluded.csv \
        && $(wc -l PrevSimulationsToBeIncluded.csv | awk ' { print $1 } ') -gt 0  ]] ; then
      for sim in $(cat PrevSimulationsToBeIncluded.csv) ; do
        for iw in $(seq 1 $(( $wave - 1 ))) ; do # Loop on waves needed to avoid taking info in WAVEN_XXXX
          grep "$sim" "WAVE$iw"/Par1D_Wave*.asc  >>  "WAVE$wave/Params.asc" || true
          grep "$sim" "WAVE$iw"/metrics_WAVE*csv >>  "WAVE$wave/Metrics.csv" || true
        done
      done
    fi
  esac

  cmd='Rscript --vanilla htune_csv2Rdata.R "'$wave'" -dir "WAVE'$wave'" $seed -par Params.asc -sim Metrics.csv'
  echo $cmd
  eval $cmd &> log/csv2Rdata.log$$
}

function history_matching() {
  echo '----------------------------------------------------------------------'
  echo "5.3 Running htune_run_$wave.sh , log: $(pwd)/log/htune_wave$wave.log$$"
  echo '----------------------------------------------------------------------'
  echo "Reading simulated metrics in WAVE${wave}/Wave${wave}_SCM.Rdata"
  echo "and corres. targets and tolerances in WAVE${wave}/Wave${wave}_REF.Rdata"

  local emulator
  emulator="htune_Emulating_Multi_Metric_Multi_LHS_new.R"  # previously "htune_Emulating_Multi_Metric_Multi_LHS.R"

  echo "source env.sh" > "htune_run_$wave.sh"
  echo "time Rscript $emulator -wave $wave -tau $tau \
     -cutoff $cutoff -sample_size $sample_size -sample_size_next_design \
     $sample_size_next_design $seed" >> "htune_run_$wave.sh"
  log_="$(pwd)/log/htune_wave$wave.log$$"
  run time bash "htune_run_$wave.sh" &> $log_
  if [ $? -ne 0 ]; then
       if [ "$(  grep 'lm.formula.* \~ [0-9]' $log )" != "" ] ; then 
           echo no parameter dependency found by htexplo for the folowing metrics
           echo You should remove them and rerun
           grep 'lm.formula.* \~ [0-9]' $log
       fi
       echo The history matching failed. Log in $log
       exit 1
  fi

  run gs -sDEVICE=pdfwrite -dPDFSETTINGS=/ebook -q -o "InputSpace_wave${wave}_light.pdf" "InputSpace_wave${wave}.pdf"
  mv "InputSpace_wave${wave}"*pdf "WAVE$wave/"
}

function compute_scores_for_next_wave() {
  echo '----------------------------------------------------------------------'
  echo "5.4 Computing scores for the next wave, log in log/scores_$wave.log$$"
  echo '----------------------------------------------------------------------'

  function comp_scores() {
    ./post_scores.sh 1 "$wave" -nograph &>> "$(pwd)/log/scores_$wave.log$$"
    cat score*[0-9].csv | sed -e /MAX/d | awk -F, ' { if ( $NF <= 2. ) { print $1 } } ' > BestSimulations0-2.csv
    cat score*[0-9].csv | sed -e /MAX/d | awk -F, ' { if ( $NF >= 2 && $NF <= 3. ) { print $1 } } ' > PrevSimulationsToBeIncluded.csv
    echo 'Simulations with err/tolerance < 2 :'
    local iw
    for iw in $(seq 1 "$wave") ; do
       echo "$(grep -c "\-$iw\-" BestSimulations0-2.csv)" "in wave $iw"
    done
  }
  run comp_scores
}

function do_graphics() {
  echo '======================================================================'
  echo "6/ Graphics"
  echo '======================================================================'

  local nmetrics nparams NROYfrac ymin
  nmetrics=$(echo "$metrics" | wc -w)
  nparams=$(wc -l "$param" | awk ' { print $1 } ')

  if python scatter_score.py &>> log/graph.log$$ ; then
    echo "New metrics scatter plot available in $(pwd)/score_metrics.pdf"
  else
    echo "Pb with scatter_score.py, log in log/graph.log$$"
  fi

  NROYfrac=$( cat "Remaining_space_after_wave_$wave.txt" )
  echo "Remaining fraction of initial NROY at wave $wave : $NROYfrac"
  ymin=$( echo "$NROYfrac" | awk ' { print 10^(int(log($1)/log(10))-1) } ' )
  python plot_NROY.py -cv "${cutoff_val[@]}" -cw "${cutoff_wav[@]}" -wn "$wave" -min "$ymin"

  function combine_pdf() {
    local nsubplots i if ii
    if [[ $nmetrics -gt 1 ]] ; then
      pdfjam --nup "${nmetrics}x${nmetrics}" "WAVE$wave/Plots_Metrics.pdf" --outfile tmp.pdf
    else
      cp -f "WAVE$wave/Plots_Metrics.pdf" tmp.pdf
    fi
    pdfjam "InputSpace_wave$wave.pdf" "WAVE$wave/Plots_LOO.pdf" tmp.pdf --landscape --outfile "WAVE$wave/synthesis.pdf"

    nsubplots=$(echo "$nmetrics" | awk ' { print int(($1-1)^0.5) + 1 } ') ; # echo $nsubplots
    rm -f tmp*pdf ; ii=1
    for i in $(seq -w 1 "$nparams") ; do
      (( if = $ii + $nmetrics - 1 ))
      pdfjam --nup "${nsubplots}x${nsubplots}" "WAVE$wave/Plots_Metrics.pdf" ${ii}-${if} --outfile "tmp$i.pdf" --landscape
      (( ii = $if + 1 ))
    done

    pdfjam tmp*.pdf --outfile "WAVE$wave/PlotMetrics.pdf" --landscape

  }

  run combine_pdf &> out$$ &

  echo '======================================================================'
  echo "Ending graphics"
  echo '======================================================================'
}

function set_wave_cutoff() {
  # Imposing a decreasing cutoff for the NROY definition (rather arbitrary)
  # values specified in cutoff_val and cutoff_wav
  for iiw_ in $( seq 1 "$ncutoff" ) ; do
    iiw=$(( $iiw_ - 1 ))
    #echo wave $wave iiw $iiw ${cutoff_wav[$iiw]}
    if [[ $wave -ge ${cutoff_wav[$iiw]} ]] ; then
      cutoff=${cutoff_val[$iiw]}
      #echo wave $wave iiw $iiw ${cutoff_wav[$iiw]} cutoff=$cutoff
    fi
  done
  echo "Cutoff for wave $wave : $cutoff"
}

# If sourced: returns, else run main script
if [[ ! "${BASH_SOURCE[0]}" = "$0" ]]; then return 0; fi

read_arguments_and_defaults "$@"

echo '====================================================================='
echo "1/ experiment setup and controls"
echo '====================================================================='

if [[ $GCM = "post" ]] ; then
  run_post
elif [[ -d $expdir/$wdir ]] ; then
  run_resume
else
  setup_new_experiment
fi

echo "STARTING LOOP ON WAVES ON $expdir/$wdir"

cd "$expdir/$wdir"
for wave in $waves ; do

  echo '======================================================================'
  echo "In loop on waves $wave"
  echo '======================================================================'

  set_wave_cutoff

  if [[ $GCM = "no" || $GCM = "pre" ]] ; then
    build_wave_design
    run_SCM_for_wave
    compute_metrics
    launch_postproc
  fi # GCM

  if [[ $GCM = "no" || $GCM = "post" ]] ; then
     echo '======================================================================'
     echo "5/ Building emulators and estimating NROY space $wave"
     echo '======================================================================'

    update_sample_size
    use_previous_waves_results
    history_matching
    compute_scores_for_next_wave
    if [[ $do_plot -eq 1 ]] ; then do_graphics; fi
  fi

  echo '======================================================================'
  echo "Ending wave $wave"
  echo '======================================================================'
done

echo "End of $0"