diff --git a/examples/infeas_analysis_cluster.ipynb b/examples/infeas_analysis_cluster.ipynb
index bb54ddd..ab2f4ea 100644
--- a/examples/infeas_analysis_cluster.ipynb
+++ b/examples/infeas_analysis_cluster.ipynb
@@ -424,7 +424,8 @@
}
],
"source": [
- "campaign, samples = analyse.read_campaign(campaign_name=\"example_cluster\")\n",
+ "campaign = analyse.get_campaign(campaign_name=\"example_cluster\")\n",
+ "samples = analyse.get_samples(campaign)\n",
"samples"
]
},
diff --git a/examples/infeas_analysis_local.ipynb b/examples/infeas_analysis_local.ipynb
index be53169..0885f2a 100644
--- a/examples/infeas_analysis_local.ipynb
+++ b/examples/infeas_analysis_local.ipynb
@@ -411,7 +411,8 @@
}
],
"source": [
- "campaign, samples = analyse.read_campaign(campaign_name=\"example_local\")\n",
+ "campaign = analyse.get_campaign(campaign_name=\"example_local\")\n",
+ "samples = analyse.get_samples(campaign)\n",
"samples"
]
},
diff --git a/infeas/analyse.py b/infeas/analyse.py
index 032923d..0817813 100644
--- a/infeas/analyse.py
+++ b/infeas/analyse.py
@@ -9,8 +9,8 @@
from infeas.eval import QOIS, WORK_DIR
-def read_campaign(campaign_name: str) -> pd.DataFrame:
- """Read in evaluated campaign and return dataframe.
+def get_campaign(campaign_name: str) -> uq.campaign.Campaign:
+ """Read in evaluated campaign and return campaign object.
Fetches the latest campaign with matching name.
@@ -23,8 +23,6 @@ def read_campaign(campaign_name: str) -> pd.DataFrame:
-------
uq.campaign.Campaign
easyVVUQ campaign
- pd.DataFrame
- Evaluation results
Raises
------
@@ -52,13 +50,27 @@ def read_campaign(campaign_name: str) -> pd.DataFrame:
db_location=db_location_prefixed, name=campaign_name, work_dir=WORK_DIR
)
+ return campaign
+
+
+def get_samples(campaign: uq.campaign.Campaign) -> pd.DataFrame:
+ """Get samples from campaign object.
+
+ Parameters
+ ----------
+ campaign : uq.campaign.Campaign
+ campaign object
+
+ Returns
+ -------
+ pd.DataFrame
+ evaluated samples
+ """
samples = campaign.get_collation_result()
- sample_count = samples.shape[0]
- print(f"Campaign read in. Number of samples = {sample_count}")
- # Drop strange multi-index of 0
+ # Drop unnecessary multi-index of 0
samples.columns = samples.columns.droplevel(1)
- return campaign, samples
+ return samples
def describe_qois(results: pd.DataFrame) -> pd.DataFrame:
diff --git a/infeas/decoder.py b/infeas/decoder.py
index c774a41..431fe1e 100644
--- a/infeas/decoder.py
+++ b/infeas/decoder.py
@@ -73,27 +73,21 @@ def _process_raw_data(self, raw_data: Dict[str, Any]) -> Dict[str, float]:
for ineq_constrs_key in ineq_constrs_keys
}
- # Only want violated constraint values
- # Coerce feasible inequality constraints (> 0) = 0.0
- # TODO Not sure if we want to mask non-violated constraint
- # values at this stage: infeasibile responses only
- vio_ineq_constrs_dict = {}
- for key, value in ineq_constrs_dict.items():
- if value > 0:
- vio_ineq_constrs_dict[key] = 0.0
- else:
- vio_ineq_constrs_dict[key] = value
-
- # Merge individual eq and ineq violated constraint values
- responses = responses | eq_constrs_dict | vio_ineq_constrs_dict
+ # Merge individual eq and ineq constraint values
+ responses = responses | eq_constrs_dict | ineq_constrs_dict
# Calculate RMS constraint residuals for violated constraints only
# Create arrays from constraints dicts
- eq_constrs = np.array(list(eq_constrs_dict.values()))
- vio_ineq_constrs = np.array(list(vio_ineq_constrs_dict.values()))
- vio_constrs = np.concatenate((eq_constrs, vio_ineq_constrs))
- rms_vio_constr_res = np.sqrt(np.mean(vio_constrs**2))
- responses["rms_vio_constr_res"] = rms_vio_constr_res
+ # Process satisfied con: c > 0
+ # Everyone else satisfied con: c < 0
+ # Flip sign to agree with standard
+ ineq_constrs = -np.array(list(ineq_constrs_dict.values()))
+
+ # TODO Would be better to override method in notebook for testing
+
+ # Use w metric: most violated inequality constraint
+ w = np.max(ineq_constrs)
+ responses["w"] = w
return responses
diff --git a/infeas/mfile_to_in.py b/infeas/mfile_to_in.py
index f5867bc..2b5da5b 100644
--- a/infeas/mfile_to_in.py
+++ b/infeas/mfile_to_in.py
@@ -9,6 +9,7 @@
from pathlib import Path
from typing import Optional, Sequence
import re
+from shutil import copy
"""Unfortunately, f-values have to be detected in two ways. This is because
@@ -120,34 +121,111 @@
def convert(
- mfile_name: str,
original_in_name: str,
- sol_in_name: str,
+ new_in_name: str,
no_optimisation: Optional[bool] = True,
n_equalities: Optional[int] = None,
+ mfile_name: Optional[str] = None,
+ remove_f_value_inits: Optional[bool] = False,
) -> Path:
- """Convert mfile to input file, preserving optimisation parameter vector.
+ """Convert input file to a modified new one.
+
+ If input file and n_equalities provided, just convert input to inequalities.
+ If mfile is provided, copy solution vector to new input file.
Parameters
----------
- mfile_name : str
- name of mfile to convert
original_in_name : str
the IN.DAT used to create the MFILE.DAT solution
- sol_in_name : str
- name of input file to be created at solution vector
+ new_in_name : str
+ modified input file to be created
no_optimisation : Optional[bool], optional
convert IN.DAT to non-optimising run, by default True
n_equalities : Optional[int], optional
how many of the constraints to make equalities. The first n_equalities
constraints will be left as equalities, the rest converted to
inequalities, by default None
+ mfile_name : Optional[str], optional
+ if mfile supplied, extract solution vector into new input file, by
+ default None
+ remove_f_value_inits : Optional[bool], optional
+ remove f-value initialisations, e.g. `fthresh = 0.8`, by default False
Returns
-------
Path
path to solution input file
"""
+ # Create new input file, optionally overwriting with solution vector
+ if mfile_name:
+ # Extract solution vector from mfile, copy into new input file
+ sol_in_dat_path = overwrite_sol_vector(
+ mfile_name=mfile_name,
+ original_in_name=original_in_name,
+ new_in_name=new_in_name,
+ )
+ else:
+ # Just copy original input to new path
+ sol_in_dat_path = copy(original_in_name, new_in_name)
+
+ # Additional modifications to top of input file
+ top_content = []
+ if no_optimisation:
+ top_content.extend(["* Once through only: no optimisation", "ioptimz = -2"])
+
+ if n_equalities is not None:
+ # Write neqns count
+ top_content.extend(
+ [
+ "* Define number of equality constraints, and",
+ "* use inequality constraints: corresponding f-values removed",
+ f"neqns = {n_equalities}",
+ ]
+ )
+
+ new_content = [line + "\n" for line in top_content]
+
+ # Write top and main content, removing f-values if necessary
+ with open(sol_in_dat_path, "r") as f:
+ content = f.readlines()
+
+ if n_equalities:
+ # Remove all f-values as opt params
+ main_content = remove_f_values(content, remove_f_value_inits)
+ else:
+ main_content = content
+
+ new_content.extend(main_content)
+
+ with open(sol_in_dat_path, "w") as f:
+ f.writelines(new_content)
+
+ return sol_in_dat_path
+
+
+def overwrite_sol_vector(
+ mfile_name: str, original_in_name: str, new_in_name: str
+) -> Path:
+ """Extract solution vector from mfile, copy into new input file.
+
+ Parameters
+ ----------
+ mfile_name : str
+ output containing solution vector
+ original_in_name : str
+ original input file
+ new_in_name : str
+ input file with initial optimisation parameter vector overwritten with
+ the solution vector
+
+ Returns
+ -------
+ Path
+ path to new solution vector input file
+ """
+ # Pattern for single array element, e.g. fimp(13)
+ array_el_re = re.compile(r"(\w+)\((\d+)\)")
+
# First, extract the solution (optimisation parameters) from the mfile
# Create Mfile object from mfile
mfile_path = Path(mfile_name)
@@ -189,63 +267,53 @@ def convert(
continue
# Otherwise update new IN.DAT value with solution value
- in_dat.add_parameter(var_name, value)
+ matches = array_el_re.match(var_name)
+ if matches is not None:
+ # Found a single array element e.g. fimp(13)
+ array_name = matches.group(1)
+ array_index = matches.group(2)
+ py_array_index = int(array_index) - 1
+
+ in_dat.change_array(
+ array_name=array_name, array_index=py_array_index, array_value=value
+ )
+ else:
+ # Set individual parameter
+ in_dat.add_parameter(var_name, value)
- print(f"Solution f-values ignored = {dropped_sol_f_value_count}")
+ print(f"Optimisation parameter f-values ignored = {dropped_sol_f_value_count}")
# Write out new IN.DAT, with optimisation parameters set to original
# solution vector
- sol_in_dat_path = Path(sol_in_name)
+ sol_in_dat_path = Path(new_in_name)
in_dat.write_in_dat(sol_in_dat_path)
-
- # Additional modifications to top of input file
- top_content = []
- if no_optimisation:
- top_content.extend(["* Once through only: no optimisation", "ioptimz = -2"])
-
- if n_equalities is not None:
- # Write neqns count
- top_content.extend(
- [
- "* Define number of equality constraints, and",
- "* use inequality constraints: corresponding f-values removed",
- f"neqns = {n_equalities}",
- ]
- )
-
- new_content = [line + "\n" for line in top_content]
-
- # Write top and main content, removing f-values if necessary
- with open(sol_in_dat_path, "r+") as f:
- content = f.readlines()
- f.seek(0, 0)
-
- if n_equalities:
- # Remove all f-values
- main_content = remove_f_values(content)
- else:
- main_content = content
-
- new_content.extend(main_content)
- f.writelines(new_content)
-
return sol_in_dat_path
-def remove_f_values(lines_with_f_values: Sequence[str]) -> list[str]:
+def remove_f_values(
+ lines_with_f_values: Sequence[str], remove_f_value_inits: bool
+) -> list[str]:
"""Remove f-value optimisation parameters from input file lines.
+ Removes `ixc = 23 * fthresh` line if opt param 23 is f-value.
+ Optionally removes f-value value initialisations `fthresh = 0.8`.
+
Parameters
----------
lines_with_f_values : Sequence[str]
lines from original input file
+ remove_f_value_inits : bool
+ remove f-value initialisations, e.g. `fthresh = 0.8`
Returns
-------
list[str]
lines with f-value optimisation parameters removed
"""
+ # f-value opt params removed (ixc = 23 * fdene)
f_value_removal_count = 0
+ # f-value initialisations removed (fthresh = 0.8)
+ f_value_init_removal_count = 0
# Lines that will be included in new IN.DAT
lines_without_f_values = []
# Opt params beginning with f that aren't on the ignore list
@@ -255,10 +323,14 @@ def remove_f_values(lines_with_f_values: Sequence[str]) -> list[str]:
# Don't also search for comment (e.g. ixc = 23 * fdene) as it may not be there:
# Just look for opt param number
opt_param_re = re.compile(r"ixc\s*=\s*(\d+)")
+
# If not already found in the f-value opt param number list, attempt to
# find any opt params starting with f (depends on "*" comment)
opt_param_beginning_with_f_re = re.compile(r"ixc\s*=\s*(\d+)\s*\*\s*(f\w+)")
+ # Any line starting with f (could be f-value value definition)
+ line_starting_with_f_re = re.compile(r"^(f\w*)")
+
for line in lines_with_f_values:
matches = opt_param_re.match(line)
if matches is not None:
@@ -272,9 +344,6 @@ def remove_f_values(lines_with_f_values: Sequence[str]) -> list[str]:
f_value_removal_count += 1
continue
- # Non f-value opt param or other line: keep
- lines_without_f_values.append(line)
-
# Check if line contains opt param beginning with an f:
# Not known f-value (would be in numbers list), but suspicious as could be
matches = opt_param_beginning_with_f_re.match(line)
@@ -285,6 +354,20 @@ def remove_f_values(lines_with_f_values: Sequence[str]) -> list[str]:
# Could be an f-value
suspect_opt_params_beginning_with_f.append(line)
+ # Optionally, remove any f-value initialisations
+ # e.g. `fthresh = 0.8`
+ if remove_f_value_inits:
+ matches = line_starting_with_f_re.match(line)
+ if matches is not None:
+ if matches.group(1) not in STARTS_WITH_F_BUT_NOT_F_VALUE_NAMES:
+ # Found f-value initialisation: drop the line
+ f_value_init_removal_count += 1
+ continue
+
+ # Non f-value opt param, (optionally not) f-value initialisation or other
+ # line: keep
+ lines_without_f_values.append(line)
+
# Print out potential f-value vars if any
if len(suspect_opt_params_beginning_with_f) > 0:
raise ValueError(
@@ -293,4 +376,7 @@ def remove_f_values(lines_with_f_values: Sequence[str]) -> list[str]:
else:
print(f"f-values removed as optimisation parameters = {f_value_removal_count}")
+ if remove_f_value_inits:
+ print(f"f-value initialisations removed = {f_value_init_removal_count}")
+
return lines_without_f_values
diff --git a/max_net_elec/max_net_elec_analysis.ipynb b/max_net_elec/max_net_elec_analysis.ipynb
index 94ee8b4..ee8f2c9 100644
--- a/max_net_elec/max_net_elec_analysis.ipynb
+++ b/max_net_elec/max_net_elec_analysis.ipynb
@@ -1247,7 +1247,8 @@
}
],
"source": [
- "campaign, results = analyse.read_campaign(campaign_name)\n",
+ "campaign = analyse.get_campaign(campaign_name)\n",
+ "results = analyse.get_samples(campaign)\n",
"results"
]
},
diff --git a/min_rmajor/min_rmajor_analysis.ipynb b/min_rmajor/min_rmajor_analysis.ipynb
index 9551910..91d4984 100644
--- a/min_rmajor/min_rmajor_analysis.ipynb
+++ b/min_rmajor/min_rmajor_analysis.ipynb
@@ -1246,7 +1246,8 @@
}
],
"source": [
- "campaign, results = analyse.read_campaign(campaign_name)\n",
+ "campaign = analyse.get_campaign(campaign_name)\n",
+ "results = analyse.get_samples(campaign)\n",
"results"
]
},
diff --git a/ra/metrics/README.md b/ra/metrics/README.md
new file mode 100644
index 0000000..25ce4da
--- /dev/null
+++ b/ra/metrics/README.md
@@ -0,0 +1,5 @@
+# Local reliability analysis: different reliability metrics
+
+Run a local reliability analysis, comparing different reliability metrics.
+
+Originally, the requirement violation severity function $s$ was defined by the RMS violated inequality constraints. This is modified to a normalised version; normalised to the value of the severity at the solution point. The severity function then becomes a test of whether the sample point is more or less violated than the solution point.
\ No newline at end of file
diff --git a/ra/metrics/lt_ineqs_after_ra_mods.template b/ra/metrics/lt_ineqs_after_ra_mods.template
new file mode 100644
index 0000000..8b0e875
--- /dev/null
+++ b/ra/metrics/lt_ineqs_after_ra_mods.template
@@ -0,0 +1,449 @@
+* Once through only: no optimisation
+ioptimz = -2
+* Define number of equality constraints, and
+* use inequality constraints: corresponding f-values removed
+neqns = 3
+
+*--------------------------------------------------*
+
+
+*---------------Constraint Equations---------------*
+
+icc = 1 * Beta
+icc = 2 * Global power balance
+icc = 11 * Radial build
+icc = 5 * Density upper limit
+icc = 8 * Neutron wall load upper limit
+icc = 9 * Fusion power upper limit
+icc = 13 * Burn time lower limit
+icc = 15 * LH power threshold limit
+icc = 30 * Injection power upper limit
+icc = 16 * Net electric power lower limit
+icc = 24 * Beta upper limit
+icc = 25 * Peak toroidal field upper limit
+icc = 26 * Central solenoid EOF current density upper limit
+icc = 27 * Central solenoid BOP current density upper limit
+icc = 33 * I_op
+icc = 34 * Dump voltage upper limit
+icc = 35 * J_winding pack
+icc = 36 * TF coil temperature margin lower limit
+icc = 60 * Central solenoid temperature margin lower limit
+icc = 62 * taup
+icc = 65 * Dump time set by VV loads
+icc = 72 * central solenoid shear stress limit
+icc = 81 * Ne
+icc = 68 * Psep
+icc = 31 * TF coil case stress upper limit
+icc = 32 * TF coil conduit stress upper limit
+
+*---------------Iteration Variables----------------*
+
+ixc = 2 * bt
+ixc = 3 * rmajor
+boundl(3) = 8.0
+boundu(3) = 9.0
+ixc = 4 * te
+boundu(4) = 100.0
+ixc = 5 * beta
+ixc = 6 * dene
+ixc = 10 * hfact
+boundu(10) = 1.2
+ixc = 13 * tfcth
+boundl(13) = 0.7
+ixc = 16 * ohcth
+boundl(16) = 0.3
+ixc = 18 * q
+boundl(18) = 3.0
+ixc = 29 * bore
+boundl(29) = 0.1
+ixc = 37 * coheof
+ixc = 41 * fcohbop
+ixc = 44 * fvsbrnni
+ixc = 56 * tdmptf
+ixc = 57 * thkcas
+ixc = 58 * thwcndut
+boundl(58) = 0.008
+ixc = 59 * fcutfsu
+boundl(59) = 0.50
+boundu(59) = 0.94
+ixc = 60 * cpttf
+boundl(60) = 65000.0
+boundu(60) = 90000.0
+ixc = 109 * ralpne
+boundu(109) = 0.1
+ixc = 122 * oh_steel_frac
+ixc = 135 * fimp(13)
+ixc = 140 * dr_tf_wp
+boundl(140) = 0.4
+
+*---------------Cs Fatigue Variables---------------*
+
+
+*------------------- Costs 1990--------------------*
+
+
+*------------------- Costs 2015--------------------*
+
+
+*-----------------Blanket Library------------------*
+
+
+*----------------------Build-----------------------*
+
+
+*-----------------Build Variables------------------*
+
+blnkith = 0.7 * inboard blanket thickness (m); (calculated if `blktmodel>0`) (=0;0 if `iblnkith=0`)
+blnkoth = 1.0 * outboard blanket thickness (m); calculated if `blktmodel>0`
+bore = 2.003843190236783 * central solenoid inboard radius (m) (`iteration variable 29`)
+ddwex = 0.15 * cryostat thickness (m)
+d_vv_in = 0.3 * vacuum vessel inboard thickness (TF coil / shield) (m)
+d_vv_out = 0.3 * vacuum vessel outboard thickness (TF coil / shield) (m)
+d_vv_top = 0.3 * vacuum vessel topside thickness (TF coil / shield) (m) (= d_vv_bot if double-null)
+d_vv_bot = 0.3 * vacuum vessel underside thickness (TF coil / shield) (m)
+gapds = 0.02 * gap between inboard vacuum vessel and thermal shield (m) (`iteration variable 61`)
+ohcth = 0.546816593988753 * Central solenoid thickness (m) (`iteration variable 16`)
+scrapli = 0.25 * Gap between plasma and first wall; inboard side (m) (if `iscrp=1`)
+scraplo = 0.25 * Gap between plasma and first wall; outboard side (m) (if `iscrp=1`)
+shldith = 0.3 * inboard shield thickness (m) (`iteration variable 93`)
+shldoth = 0.800 * outboard shield thickness (m) (`iteration variable 94`)
+tfcth = 1.2 * inboard TF coil thickness; (centrepost for ST) (m)
+thshield_ib = 0.050 * TF-VV thermal shield thickness; inboard (m)
+thshield_ob = 0.050 * TF-VV thermal shield thickness; outboard (m)
+thshield_vb = 0.050 * TF-VV thermal shield thickness; vertical build (m)
+vvblgap = 0.02 * gap between vacuum vessel and blanket (m)
+
+*---------------Buildings Variables----------------*
+
+
+*-----------------Ccfe Hcpb Module-----------------*
+
+
+*---------------Const And Precisions---------------*
+
+
+*--------------------Constants---------------------*
+
+
+*---------------Constraint Variables---------------*
+
+bmxlim = 14.0 * maximum peak toroidal field (T) (`constraint equation 25`)
+fbetatry = 0.5 * f-value for beta limit (`constraint equation 24`; `iteration variable 36`)
+fdene = 1.2 * f-value for density limit (`constraint equation 5`; `iteration variable 9`)
+fiooic = 0.65 * f-value for TF coil operating current / critical current ratio
+fjohc = 0.6 * f-value for central solenoid current at end-of-flattop
+fjohc0 = 0.6 * f-value for central solenoid current at beginning of pulse
+fjprot = 1.0 * f-value for TF coil winding pack current density
+foh_stress = 1.0 * f-value for Tresca yield criterion in Central Solenoid
+fmaxvvstress = 1.0 * f-value for maximum permitted stress of the VV
+fvdump = 1.0 * f-value for dump voltage (`constraint equation 34`; `iteration variable 51`)
+fwalld = 1.0 * f-value for maximum wall load (`constraint equation 8`; `iteration variable 14`)
+pnetelin = 400.0 * required net electric power (MW) (`constraint equation 16`)
+powfmax = 3000 * maximum fusion power (MW) (`constraint equation 9`)
+psepbqarmax = 10.0 * maximum ratio of Psep*Bt/qAR (MWT/m) (`constraint equation 68`)
+tbrnmn = 7200.0 * minimum burn time (s) (KE - no longer itv;; see issue #706)
+walalw = 2.0 * allowable neutron wall-load (MW/m2) (`constraint equation 8`)
+taulimit = 5.0 * Lower limit on taup/taueff the ratio of alpha particle to energy confinement
+
+*-------------------Constraints--------------------*
+
+
+*------------------Cost Variables------------------*
+
+cfactr = 0.80 * Total plant availability fraction; input if `iavail=0`
+cost_model = 0 * Switch for cost model;
+iavail = 0 * Switch for plant availability model;
+output_costs = 1 * Switch for costs output;
+
+*----------------------Costs-----------------------*
+
+
+*-------------Current Drive Variables--------------*
+
+bscfmax = 0.95 * maximum fraction of plasma current from bootstrap; if `bscfmax < 0`;
+etaech = 0.5 * ECH wall plug to injector efficiency
+gamma_ecrh = 0.30 * User input ECRH gamma (1;0e20 A/(W m^2))
+iefrf = 10 * Switch for current drive efficiency model;
+pheat = 75.0 * heating power not used for current drive (MW) (`iteration variable 11`)
+pinjalw = 200.0 * maximum allowable value for injected power (MW) (`constraint equation 30`)
+
+*-------------------Dcll Module--------------------*
+
+
+*------------Define Iteration Variables------------*
+
+
+*----------------Divertor Variables----------------*
+
+divfix = 0.62 * divertor structure vertical thickness (m)
+
+*------------------Error Handling------------------*
+
+
+*-------------------Final Module-------------------*
+
+
+*-------------------Fson Library-------------------*
+
+
+*-------------------Fson Path M--------------------*
+
+
+*------------------Fson String M-------------------*
+
+
+*-------------------Fson Value M-------------------*
+
+
+*----------------Function Evaluator----------------*
+
+
+*--------------------Fw Module---------------------*
+
+
+*-------------------Fwbs Module--------------------*
+
+
+*------------------Fwbs Variables------------------*
+
+inuclear = 1 * switch for nuclear heating in the coils;
+qnuc = 1.3e4 * nuclear heating in the coils (W) (`inuclear=1`)
+primary_pumping = 3 * Switch for pumping power for primary coolant (mechanical power only and peak first wall
+secondary_cycle = 2 * Switch for power conversion cycle;
+vfshld = 0.60 * coolant void fraction in shield
+etaiso = 0.9 * isentropic efficiency of FW and blanket coolant pumps
+etahtp = 0.87 * electrical efficiency of primary coolant pumps
+
+*-----------------Global Variables-----------------*
+
+runtitle = generic large tokamak * short descriptive title for the run
+
+*-------------Heat Transport Variables-------------*
+
+etath = 0.4 * thermal to electric conversion efficiency if `secondary_cycle=2`; otherwise calculated;
+ipowerflow = 0 * switch for power flow model;
+iprimshld = 1 * Switch for shield thermal power destiny;
+
+*--------------------Ife Module--------------------*
+
+
+*------------------Ife Variables-------------------*
+
+
+*------------Impurity Radiation Module-------------*
+
+coreradius = 0.75 * coreradius /0;6/ ; normalised radius defining the 'core' region
+coreradiationfraction = 0.6 * coreradiationfraction /1;0/ ; fraction of radiation from 'core' region that is subtracted from the loss power
+fimp(1) = 0.9
+fimp(2) = 0.1
+fimp(3) = 0.0
+fimp(4) = 0.0
+fimp(5) = 0.0
+fimp(6) = 0.0
+fimp(7) = 0.0
+fimp(8) = 0.0
+fimp(9) = 0.0
+fimp(10) = 0.0
+fimp(11) = 0.0
+fimp(12) = 0.0
+fimp(13) = 0.000597755323387789
+fimp(14) = 5e-06
+
+*-------------------Init Module--------------------*
+
+
+*-------------------Main Module--------------------*
+
+
+*------------------Maths Library-------------------*
+
+
+*--------------Neoclassics Constants---------------*
+
+
+*----------------Neoclassics Module----------------*
+
+
+*---------------------Numerics---------------------*
+
+minmax = 1 *
+neqns = 3 * neqns /0/ ; number of equality constraints to be satisfied
+epsvmc = 1e-7 * epsvmc /1;0e-6/ ; error tolerance for VMCON
+
+*----------------Pf Power Variables----------------*
+
+
+*------------------Pfcoil Module-------------------*
+
+
+*-----------------Pfcoil Variables-----------------*
+
+alstroh = 7.5d8 * allowable hoop stress in Central Solenoid structural material (Pa)
+coheof = 21443595.371072624 * Central solenoid overall current density at end of flat-top (A/m2) (`iteration variable 37`) (`sweep variable 62`)
+cptdin = 4.0d4, 4.0d4, 4.0d4, 4.0d4, 4.0d4, 4.0d4, 4.0d4, 4.0d4 * peak current per turn input for PF coil i (A)
+fcohbop = 0.93491189654662 * ratio of central solenoid overall current density at beginning of pulse / end of flat-top
+fcuohsu = 0.70 * copper fraction of strand in central solenoid
+ipfloc = 2,2,3,3 * Switch for location of PF coil group i;
+isumatoh = 1 * switch for superconductor material in central solenoid;
+isumatpf = 3 * switch for superconductor material in PF coils;
+ncls = 1,1,2,2 * number of PF coils in group j
+ngrp = 4 * number of groups of PF coils; Symmetric coil pairs should all be in the same group
+ohhghf = 0.9 * Central solenoid height / TF coil internal height
+oh_steel_frac = 0.4856940627014451 * central solenoid steel fraction (`iteration variable 122`)
+rjconpf = 1.1d7, 1.1d7, 6.d6, 6.d6, 8.d6, 8.0d6, 8.0d6, 8.0d6 * average winding pack current density of PF coil i (A/m2) at time of peak
+rpf2 = -1.825 * offset (m) of radial position of `ipfloc=2` PF coils from being at
+sigpfcf = 0.666 * fraction of JxB hoop force supported by steel case for superconducting PF coils (`ipfres=0`)
+zref(1) = 3.6
+zref(2) = 1.2
+zref(3) = 1.0
+zref(4) = 2.8
+zref(5) = 1.0
+zref(6) = 1.0
+zref(7) = 1.0
+zref(8) = 1.0
+zref(9) = 1.0
+zref(10) = 1.0
+
+*-------------Physics Functions Module-------------*
+
+
+*------------------Physics Module------------------*
+
+
+*----------------Physics Variables-----------------*
+
+alphan = 1.00 * density profile index
+alphat = 1.45 * temperature profile index
+aspect = 3.0 * aspect ratio (`iteration variable 1`)
+beta = 0.03230408815355488 * total plasma beta (`iteration variable 5`) (calculated if stellarator)
+bt = 5.318322174644904 * toroidal field on axis (T) (`iteration variable 2`)
+dene = 7.796223900029837e+19 * electron density (/m3) (`iteration variable 6`)
+dnbeta = 3.0 * Troyon-like coefficient for beta scaling calculated
+fgwsep = 0.5 * fraction of Greenwald density to set as separatrix density; If `<0`; separatrix
+fkzohm = 1.02 * Zohm elongation scaling adjustment factor (`ishape=2; 3`)
+fvsbrnni = 0.4242184436680697 * fraction of the plasma current produced by non-inductive means (`iteration variable 44`)
+gamma = 0.3 * Ejima coefficient for resistive startup V-s formula
+hfact = 1.185971818905028 * H factor on energy confinement times; radiation corrected (`iteration variable 10`);
+ibss = 4 * switch for bootstrap current scaling
+iculbl = 1 * switch for beta limit scaling (`constraint equation 24`)
+icurr = 4 * switch for plasma current scaling to use
+idensl = 7 * switch for density limit to enforce (`constraint equation 5`)
+ifalphap = 1 * switch for fast alpha pressure calculation
+iinvqd = 1 * switch for inverse quadrature in L-mode scaling laws 5 and 9;
+ipedestal = 1 * switch for pedestal profiles;
+neped = 0.5e20 * electron density of pedestal [m-3] (`ipedestal==1)
+nesep = 0.2e20 * electron density at separatrix [m-3] (`ipedestal==1)
+plasma_res_factor = 0.7 * plasma resistivity pre-factor
+rhopedn = 0.94 * r/a of density pedestal (`ipedestal==1`)
+rhopedt = 0.94 * r/a of temperature pedestal (`ipedestal==1`)
+tbeta = 2.0 * temperature profile index beta (`ipedestal==1)
+teped = 5.5 * electron temperature of pedestal (keV) (`ipedestal==1; ieped=0; calculated for ieped=1`)
+tesep = 0.1 * electron temperature at separatrix (keV) (`ipedestal==1`) calculated if reinke
+iprofile = 1 * switch for current profile consistency;
+isc = 34 * switch for energy confinement time scaling law (see description in `tauscl`)
+ishape = 0 * switch for plasma cross-sectional shape calculation;
+kappa = 1.85 * plasma separatrix elongation (calculated if `ishape = 1-5; 7 or 9-10`)
+q = 3.7339078193128556 * Safety factor 'near' plasma edge (`iteration variable 18`) equal to q95
+q0 = 1.0 * safety factor on axis
+ralpne = 0.060238763988650204 * thermal alpha density/electron density (`iteration variable 109`)
+rmajor = 8.0 * plasma major radius (m) (`iteration variable 3`)
+i_single_null = 1 * switch for single null / double null plasma;
+ssync = 0.6 * synchrotron wall reflectivity factor
+te = 12.221383528378944 * volume averaged electron temperature (keV) (`iteration variable 4`)
+triang = 0.5 * plasma separatrix triangularity (calculated if `ishape = 1; 3-5 or 7`)
+
+*----------------------Power-----------------------*
+
+
+*------------Primary Pumping Variables-------------*
+
+
+*------------------Process Input-------------------*
+
+
+*------------------Process Output------------------*
+
+
+*-----------------Profiles Module------------------*
+
+
+*-----------------Pulse Variables------------------*
+
+lpulse = 1 * Switch for reactor model;
+
+*-----------------Rebco Variables------------------*
+
+
+*------------------Reinke Module-------------------*
+
+
+*-----------------Reinke Variables-----------------*
+
+
+*---------------Resistive Materials----------------*
+
+
+*-------------------Scan Module--------------------*
+
+
+*-----------------Sctfcoil Module------------------*
+
+
+*----------------Startup Variables-----------------*
+
+
+*------------Stellarator Configuration-------------*
+
+
+*----------------Stellarator Module----------------*
+
+
+*--------------Stellarator Variables---------------*
+
+
+*---------------Structure Variables----------------*
+
+
+*-----------------Tfcoil Variables-----------------*
+
+sig_tf_case_max = 7.5e8 * Allowable maximum shear stress (Tresca criterion) in TF coil case (Pa)
+sig_tf_wp_max = 7.5e8 * Allowable maximum shear stress (Tresca criterion) in TF coil conduit (Pa)
+casthi = 0.06 * inboard TF coil case plasma side thickness (m) (calculated for stellarators)
+casths = 0.05 * inboard TF coil sidewall case thickness (m) (calculated for stellarators)
+cpttf = 85462.67500253802 * TF coil current per turn (A); (calculated for stellarators) (calculated for
+dhecoil = 0.01 * diameter of central helium channel in TF winding (m)
+fcutfsu = 0.8231999768826475 * copper fraction of cable conductor (TF coils)
+i_tf_sc_mat = 1 * Switch for superconductor material in TF coils;
+ripmax = 0.6 * aximum allowable toroidal field ripple amplitude at plasma edge (%)
+tdmptf = 17.97282589344206 * fast discharge time for TF coil in event of quench (s) (`iteration variable 56`)
+n_tf = 16 * Number of TF coils (default = 50 for stellarators); Number of TF coils outer legs for ST
+tftmp = 4.75 * peak helium coolant temperature in TF coils and PF coils (K)
+thkcas = 0.2816873221155309 * inboard TF coil case outer (non-plasma side) thickness (m) (`iteration variable 57`)
+dr_tf_wp = 0.5153787768966674 * radial thickness of winding pack (m) (`iteration variable 140`) (issue #514)
+thwcndut = 0.008012110032981922 * TF coil conduit case thickness (m) (`iteration variable 58`)
+tinstf = 0.008 * Thickness of the ground insulation layer surrounding (m)
+tmargmin_cs = 1.5 * minimum allowable temperature margin ; CS (K)
+tmargmin = 1.5 * minimum allowable temperature margin ; TFC AND CS (K)
+vdalw = 10.0 * max voltage across TF coil during quench (kV) (`iteration variable 52`)
+vftf = 0.3 * coolant fraction of TFC 'cable' (`i_tf_sup=1`); or of TFC leg (`i_tf_ssup=0`)
+
+*-----------------Times Variables------------------*
+
+pulsetimings = 0 * Switch for pulse timings (if lpulse=1);
+tdwell = 1800.0 * time between pulses in a pulsed reactor (s) (`iteration variable 17`)
+tramp = 500.0 * initial PF coil charge time (s); if pulsed; = tohs
+
+*--------------------Utilities---------------------*
+
+
+*-----------------Vacuum Variables-----------------*
+
+
+*--------------Water Usage Variables---------------*
+
+******************************
+* Modifications for reliability analysis: overwrites any previous values
+******************************
+fimp(14) = $fimp_14 * W impurity fraction
+psepbqarmax = $psepbqarmax * divertor protection limit
+triang = $triang * triangularity
\ No newline at end of file
diff --git a/ra/metrics/lt_ineqs_before_conv_crit_mod.template b/ra/metrics/lt_ineqs_before_conv_crit_mod.template
new file mode 100644
index 0000000..22c24fe
--- /dev/null
+++ b/ra/metrics/lt_ineqs_before_conv_crit_mod.template
@@ -0,0 +1,449 @@
+* Once through only: no optimisation
+ioptimz = -2
+* Define number of equality constraints, and
+* use inequality constraints: corresponding f-values removed
+neqns = 3
+
+*--------------------------------------------------*
+
+
+*---------------Constraint Equations---------------*
+
+icc = 1 * Beta
+icc = 2 * Global power balance
+icc = 11 * Radial build
+icc = 5 * Density upper limit
+icc = 8 * Neutron wall load upper limit
+icc = 9 * Fusion power upper limit
+icc = 13 * Burn time lower limit
+icc = 15 * LH power threshold limit
+icc = 30 * Injection power upper limit
+icc = 16 * Net electric power lower limit
+icc = 24 * Beta upper limit
+icc = 25 * Peak toroidal field upper limit
+icc = 26 * Central solenoid EOF current density upper limit
+icc = 27 * Central solenoid BOP current density upper limit
+icc = 33 * I_op
+icc = 34 * Dump voltage upper limit
+icc = 35 * J_winding pack
+icc = 36 * TF coil temperature margin lower limit
+icc = 60 * Central solenoid temperature margin lower limit
+icc = 62 * taup
+icc = 65 * Dump time set by VV loads
+icc = 72 * central solenoid shear stress limit
+icc = 81 * Ne
+icc = 68 * Psep
+icc = 31 * TF coil case stress upper limit
+icc = 32 * TF coil conduit stress upper limit
+
+*---------------Iteration Variables----------------*
+
+ixc = 2 * bt
+ixc = 3 * rmajor
+boundl(3) = 8.0
+boundu(3) = 9.0
+ixc = 4 * te
+boundu(4) = 100.0
+ixc = 5 * beta
+ixc = 6 * dene
+ixc = 10 * hfact
+boundu(10) = 1.2
+ixc = 13 * tfcth
+boundl(13) = 0.7
+ixc = 16 * ohcth
+boundl(16) = 0.3
+ixc = 18 * q
+boundl(18) = 3.0
+ixc = 29 * bore
+boundl(29) = 0.1
+ixc = 37 * coheof
+ixc = 41 * fcohbop
+ixc = 44 * fvsbrnni
+ixc = 56 * tdmptf
+ixc = 57 * thkcas
+ixc = 58 * thwcndut
+boundl(58) = 0.008
+ixc = 59 * fcutfsu
+boundl(59) = 0.50
+boundu(59) = 0.94
+ixc = 60 * cpttf
+boundl(60) = 65000.0
+boundu(60) = 90000.0
+ixc = 109 * ralpne
+boundu(109) = 0.1
+ixc = 122 * oh_steel_frac
+ixc = 135 * fimp(13)
+ixc = 140 * dr_tf_wp
+boundl(140) = 0.4
+
+*---------------Cs Fatigue Variables---------------*
+
+
+*------------------- Costs 1990--------------------*
+
+
+*------------------- Costs 2015--------------------*
+
+
+*-----------------Blanket Library------------------*
+
+
+*----------------------Build-----------------------*
+
+
+*-----------------Build Variables------------------*
+
+blnkith = 0.7 * inboard blanket thickness (m); (calculated if `blktmodel>0`) (=0;0 if `iblnkith=0`)
+blnkoth = 1.0 * outboard blanket thickness (m); calculated if `blktmodel>0`
+bore = 1.9538096403866128 * central solenoid inboard radius (m) (`iteration variable 29`)
+ddwex = 0.15 * cryostat thickness (m)
+d_vv_in = 0.3 * vacuum vessel inboard thickness (TF coil / shield) (m)
+d_vv_out = 0.3 * vacuum vessel outboard thickness (TF coil / shield) (m)
+d_vv_top = 0.3 * vacuum vessel topside thickness (TF coil / shield) (m) (= d_vv_bot if double-null)
+d_vv_bot = 0.3 * vacuum vessel underside thickness (TF coil / shield) (m)
+gapds = 0.02 * gap between inboard vacuum vessel and thermal shield (m) (`iteration variable 61`)
+ohcth = 0.5699418715863596 * Central solenoid thickness (m) (`iteration variable 16`)
+scrapli = 0.25 * Gap between plasma and first wall; inboard side (m) (if `iscrp=1`)
+scraplo = 0.25 * Gap between plasma and first wall; outboard side (m) (if `iscrp=1`)
+shldith = 0.3 * inboard shield thickness (m) (`iteration variable 93`)
+shldoth = 0.800 * outboard shield thickness (m) (`iteration variable 94`)
+tfcth = 1.2 * inboard TF coil thickness; (centrepost for ST) (m)
+thshield_ib = 0.050 * TF-VV thermal shield thickness; inboard (m)
+thshield_ob = 0.050 * TF-VV thermal shield thickness; outboard (m)
+thshield_vb = 0.050 * TF-VV thermal shield thickness; vertical build (m)
+vvblgap = 0.02 * gap between vacuum vessel and blanket (m)
+
+*---------------Buildings Variables----------------*
+
+
+*-----------------Ccfe Hcpb Module-----------------*
+
+
+*---------------Const And Precisions---------------*
+
+
+*--------------------Constants---------------------*
+
+
+*---------------Constraint Variables---------------*
+
+bmxlim = 14.0 * maximum peak toroidal field (T) (`constraint equation 25`)
+fbetatry = 0.5 * f-value for beta limit (`constraint equation 24`; `iteration variable 36`)
+fdene = 1.2 * f-value for density limit (`constraint equation 5`; `iteration variable 9`)
+fiooic = 0.65 * f-value for TF coil operating current / critical current ratio
+fjohc = 0.6 * f-value for central solenoid current at end-of-flattop
+fjohc0 = 0.6 * f-value for central solenoid current at beginning of pulse
+fjprot = 1.0 * f-value for TF coil winding pack current density
+foh_stress = 1.0 * f-value for Tresca yield criterion in Central Solenoid
+fmaxvvstress = 1.0 * f-value for maximum permitted stress of the VV
+fvdump = 1.0 * f-value for dump voltage (`constraint equation 34`; `iteration variable 51`)
+fwalld = 1.0 * f-value for maximum wall load (`constraint equation 8`; `iteration variable 14`)
+pnetelin = 400.0 * required net electric power (MW) (`constraint equation 16`)
+powfmax = 3000 * maximum fusion power (MW) (`constraint equation 9`)
+psepbqarmax = 10.0 * maximum ratio of Psep*Bt/qAR (MWT/m) (`constraint equation 68`)
+tbrnmn = 7200.0 * minimum burn time (s) (KE - no longer itv;; see issue #706)
+walalw = 2.0 * allowable neutron wall-load (MW/m2) (`constraint equation 8`)
+taulimit = 5.0 * Lower limit on taup/taueff the ratio of alpha particle to energy confinement
+
+*-------------------Constraints--------------------*
+
+
+*------------------Cost Variables------------------*
+
+cfactr = 0.80 * Total plant availability fraction; input if `iavail=0`
+cost_model = 0 * Switch for cost model;
+iavail = 0 * Switch for plant availability model;
+output_costs = 1 * Switch for costs output;
+
+*----------------------Costs-----------------------*
+
+
+*-------------Current Drive Variables--------------*
+
+bscfmax = 0.95 * maximum fraction of plasma current from bootstrap; if `bscfmax < 0`;
+etaech = 0.5 * ECH wall plug to injector efficiency
+gamma_ecrh = 0.30 * User input ECRH gamma (1;0e20 A/(W m^2))
+iefrf = 10 * Switch for current drive efficiency model;
+pheat = 75.0 * heating power not used for current drive (MW) (`iteration variable 11`)
+pinjalw = 200.0 * maximum allowable value for injected power (MW) (`constraint equation 30`)
+
+*-------------------Dcll Module--------------------*
+
+
+*------------Define Iteration Variables------------*
+
+
+*----------------Divertor Variables----------------*
+
+divfix = 0.62 * divertor structure vertical thickness (m)
+
+*------------------Error Handling------------------*
+
+
+*-------------------Final Module-------------------*
+
+
+*-------------------Fson Library-------------------*
+
+
+*-------------------Fson Path M--------------------*
+
+
+*------------------Fson String M-------------------*
+
+
+*-------------------Fson Value M-------------------*
+
+
+*----------------Function Evaluator----------------*
+
+
+*--------------------Fw Module---------------------*
+
+
+*-------------------Fwbs Module--------------------*
+
+
+*------------------Fwbs Variables------------------*
+
+inuclear = 1 * switch for nuclear heating in the coils;
+qnuc = 1.3e4 * nuclear heating in the coils (W) (`inuclear=1`)
+primary_pumping = 3 * Switch for pumping power for primary coolant (mechanical power only and peak first wall
+secondary_cycle = 2 * Switch for power conversion cycle;
+vfshld = 0.60 * coolant void fraction in shield
+etaiso = 0.9 * isentropic efficiency of FW and blanket coolant pumps
+etahtp = 0.87 * electrical efficiency of primary coolant pumps
+
+*-----------------Global Variables-----------------*
+
+runtitle = generic large tokamak * short descriptive title for the run
+
+*-------------Heat Transport Variables-------------*
+
+etath = 0.4 * thermal to electric conversion efficiency if `secondary_cycle=2`; otherwise calculated;
+ipowerflow = 0 * switch for power flow model;
+iprimshld = 1 * Switch for shield thermal power destiny;
+
+*--------------------Ife Module--------------------*
+
+
+*------------------Ife Variables-------------------*
+
+
+*------------Impurity Radiation Module-------------*
+
+coreradius = 0.75 * coreradius /0;6/ ; normalised radius defining the 'core' region
+coreradiationfraction = 0.6 * coreradiationfraction /1;0/ ; fraction of radiation from 'core' region that is subtracted from the loss power
+fimp(1) = 0.9
+fimp(2) = 0.1
+fimp(3) = 0.0
+fimp(4) = 0.0
+fimp(5) = 0.0
+fimp(6) = 0.0
+fimp(7) = 0.0
+fimp(8) = 0.0
+fimp(9) = 0.0
+fimp(10) = 0.0
+fimp(11) = 0.0
+fimp(12) = 0.0
+fimp(13) = 0.0005676083884529063
+fimp(14) = 5e-06
+
+*-------------------Init Module--------------------*
+
+
+*-------------------Main Module--------------------*
+
+
+*------------------Maths Library-------------------*
+
+
+*--------------Neoclassics Constants---------------*
+
+
+*----------------Neoclassics Module----------------*
+
+
+*---------------------Numerics---------------------*
+
+minmax = 1 *
+neqns = 3 * neqns /0/ ; number of equality constraints to be satisfied
+epsvmc = 1e-7 * epsvmc /1;0e-6/ ; error tolerance for VMCON
+
+*----------------Pf Power Variables----------------*
+
+
+*------------------Pfcoil Module-------------------*
+
+
+*-----------------Pfcoil Variables-----------------*
+
+alstroh = 7.5d8 * allowable hoop stress in Central Solenoid structural material (Pa)
+coheof = 20910063.408582743 * Central solenoid overall current density at end of flat-top (A/m2) (`iteration variable 37`) (`sweep variable 62`)
+cptdin = 4.0d4, 4.0d4, 4.0d4, 4.0d4, 4.0d4, 4.0d4, 4.0d4, 4.0d4 * peak current per turn input for PF coil i (A)
+fcohbop = 0.9253870649759701 * ratio of central solenoid overall current density at beginning of pulse / end of flat-top
+fcuohsu = 0.70 * copper fraction of strand in central solenoid
+ipfloc = 2,2,3,3 * Switch for location of PF coil group i;
+isumatoh = 1 * switch for superconductor material in central solenoid;
+isumatpf = 3 * switch for superconductor material in PF coils;
+ncls = 1,1,2,2 * number of PF coils in group j
+ngrp = 4 * number of groups of PF coils; Symmetric coil pairs should all be in the same group
+ohhghf = 0.9 * Central solenoid height / TF coil internal height
+oh_steel_frac = 0.4666244422630157 * central solenoid steel fraction (`iteration variable 122`)
+rjconpf = 1.1d7, 1.1d7, 6.d6, 6.d6, 8.d6, 8.0d6, 8.0d6, 8.0d6 * average winding pack current density of PF coil i (A/m2) at time of peak
+rpf2 = -1.825 * offset (m) of radial position of `ipfloc=2` PF coils from being at
+sigpfcf = 0.666 * fraction of JxB hoop force supported by steel case for superconducting PF coils (`ipfres=0`)
+zref(1) = 3.6
+zref(2) = 1.2
+zref(3) = 1.0
+zref(4) = 2.8
+zref(5) = 1.0
+zref(6) = 1.0
+zref(7) = 1.0
+zref(8) = 1.0
+zref(9) = 1.0
+zref(10) = 1.0
+
+*-------------Physics Functions Module-------------*
+
+
+*------------------Physics Module------------------*
+
+
+*----------------Physics Variables-----------------*
+
+alphan = 1.00 * density profile index
+alphat = 1.45 * temperature profile index
+aspect = 3.0 * aspect ratio (`iteration variable 1`)
+beta = 0.032875409159602936 * total plasma beta (`iteration variable 5`) (calculated if stellarator)
+bt = 5.3432605090882115 * toroidal field on axis (T) (`iteration variable 2`)
+dene = 7.968944559592192e+19 * electron density (/m3) (`iteration variable 6`)
+dnbeta = 3.0 * Troyon-like coefficient for beta scaling calculated
+fgwsep = 0.5 * fraction of Greenwald density to set as separatrix density; If `<0`; separatrix
+fkzohm = 1.02 * Zohm elongation scaling adjustment factor (`ishape=2; 3`)
+fvsbrnni = 0.43455894250944155 * fraction of the plasma current produced by non-inductive means (`iteration variable 44`)
+gamma = 0.3 * Ejima coefficient for resistive startup V-s formula
+hfact = 1.2 * H factor on energy confinement times; radiation corrected (`iteration variable 10`);
+ibss = 4 * switch for bootstrap current scaling
+iculbl = 1 * switch for beta limit scaling (`constraint equation 24`)
+icurr = 4 * switch for plasma current scaling to use
+idensl = 7 * switch for density limit to enforce (`constraint equation 5`)
+ifalphap = 1 * switch for fast alpha pressure calculation
+iinvqd = 1 * switch for inverse quadrature in L-mode scaling laws 5 and 9;
+ipedestal = 1 * switch for pedestal profiles;
+neped = 0.5e20 * electron density of pedestal [m-3] (`ipedestal==1)
+nesep = 0.2e20 * electron density at separatrix [m-3] (`ipedestal==1)
+plasma_res_factor = 0.7 * plasma resistivity pre-factor
+rhopedn = 0.94 * r/a of density pedestal (`ipedestal==1`)
+rhopedt = 0.94 * r/a of temperature pedestal (`ipedestal==1`)
+tbeta = 2.0 * temperature profile index beta (`ipedestal==1)
+teped = 5.5 * electron temperature of pedestal (keV) (`ipedestal==1; ieped=0; calculated for ieped=1`)
+tesep = 0.1 * electron temperature at separatrix (keV) (`ipedestal==1`) calculated if reinke
+iprofile = 1 * switch for current profile consistency;
+isc = 34 * switch for energy confinement time scaling law (see description in `tauscl`)
+ishape = 0 * switch for plasma cross-sectional shape calculation;
+kappa = 1.85 * plasma separatrix elongation (calculated if `ishape = 1-5; 7 or 9-10`)
+q = 3.6740980326674357 * Safety factor 'near' plasma edge (`iteration variable 18`) equal to q95
+q0 = 1.0 * safety factor on axis
+ralpne = 0.08184572822808084 * thermal alpha density/electron density (`iteration variable 109`)
+rmajor = 8.0 * plasma major radius (m) (`iteration variable 3`)
+i_single_null = 1 * switch for single null / double null plasma;
+ssync = 0.6 * synchrotron wall reflectivity factor
+te = 12.529208872441531 * volume averaged electron temperature (keV) (`iteration variable 4`)
+triang = 0.5 * plasma separatrix triangularity (calculated if `ishape = 1; 3-5 or 7`)
+
+*----------------------Power-----------------------*
+
+
+*------------Primary Pumping Variables-------------*
+
+
+*------------------Process Input-------------------*
+
+
+*------------------Process Output------------------*
+
+
+*-----------------Profiles Module------------------*
+
+
+*-----------------Pulse Variables------------------*
+
+lpulse = 1 * Switch for reactor model;
+
+*-----------------Rebco Variables------------------*
+
+
+*------------------Reinke Module-------------------*
+
+
+*-----------------Reinke Variables-----------------*
+
+
+*---------------Resistive Materials----------------*
+
+
+*-------------------Scan Module--------------------*
+
+
+*-----------------Sctfcoil Module------------------*
+
+
+*----------------Startup Variables-----------------*
+
+
+*------------Stellarator Configuration-------------*
+
+
+*----------------Stellarator Module----------------*
+
+
+*--------------Stellarator Variables---------------*
+
+
+*---------------Structure Variables----------------*
+
+
+*-----------------Tfcoil Variables-----------------*
+
+sig_tf_case_max = 7.5e8 * Allowable maximum shear stress (Tresca criterion) in TF coil case (Pa)
+sig_tf_wp_max = 7.5e8 * Allowable maximum shear stress (Tresca criterion) in TF coil conduit (Pa)
+casthi = 0.06 * inboard TF coil case plasma side thickness (m) (calculated for stellarators)
+casths = 0.05 * inboard TF coil sidewall case thickness (m) (calculated for stellarators)
+cpttf = 79653.5996379021 * TF coil current per turn (A); (calculated for stellarators) (calculated for
+dhecoil = 0.01 * diameter of central helium channel in TF winding (m)
+fcutfsu = 0.8279235435984497 * copper fraction of cable conductor (TF coils)
+i_tf_sc_mat = 1 * Switch for superconductor material in TF coils;
+ripmax = 0.6 * aximum allowable toroidal field ripple amplitude at plasma edge (%)
+tdmptf = 19.516156838432067 * fast discharge time for TF coil in event of quench (s) (`iteration variable 56`)
+n_tf = 16 * Number of TF coils (default = 50 for stellarators); Number of TF coils outer legs for ST
+tftmp = 4.75 * peak helium coolant temperature in TF coils and PF coils (K)
+thkcas = 0.27817131748284374 * inboard TF coil case outer (non-plasma side) thickness (m) (`iteration variable 57`)
+dr_tf_wp = 0.5446354488716441 * radial thickness of winding pack (m) (`iteration variable 140`) (issue #514)
+thwcndut = 0.008 * TF coil conduit case thickness (m) (`iteration variable 58`)
+tinstf = 0.008 * Thickness of the ground insulation layer surrounding (m)
+tmargmin_cs = 1.5 * minimum allowable temperature margin ; CS (K)
+tmargmin = 1.5 * minimum allowable temperature margin ; TFC AND CS (K)
+vdalw = 10.0 * max voltage across TF coil during quench (kV) (`iteration variable 52`)
+vftf = 0.3 * coolant fraction of TFC 'cable' (`i_tf_sup=1`); or of TFC leg (`i_tf_ssup=0`)
+
+*-----------------Times Variables------------------*
+
+pulsetimings = 0 * Switch for pulse timings (if lpulse=1);
+tdwell = 1800.0 * time between pulses in a pulsed reactor (s) (`iteration variable 17`)
+tramp = 500.0 * initial PF coil charge time (s); if pulsed; = tohs
+
+*--------------------Utilities---------------------*
+
+
+*-----------------Vacuum Variables-----------------*
+
+
+*--------------Water Usage Variables---------------*
+
+******************************
+* Modifications for reliability analysis: overwrites any previous values
+******************************
+fimp(14) = $fimp_14 * W impurity fraction
+psepbqarmax = $psepbqarmax * divertor protection limit
+triang = $triang * triangularity
\ No newline at end of file
diff --git a/ra/metrics/ra_analysis.ipynb b/ra/metrics/ra_analysis.ipynb
new file mode 100644
index 0000000..870a3b2
--- /dev/null
+++ b/ra/metrics/ra_analysis.ipynb
@@ -0,0 +1,384 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Local reliability with w severity function: analyse\n",
+ "\n",
+ "Evaluate local reliability for 3 uncertain parameters (fimp_14, psepbqarmax and triang). This uses the $w$ worst-violated inequality constraint severity metric."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 103,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "import seaborn as sns\n",
+ "import numpy as np\n",
+ "import pandas as pd\n",
+ "from infeas import analyse"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 104,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Reading in campaign database.\n",
+ "Reading in campaign database.\n"
+ ]
+ }
+ ],
+ "source": [
+ "campaign_names = [\n",
+ " \"lt_3p_before_conv_crit_mod\", # with ra mods, but without conv crit mod\n",
+ " \"lt_3p_after_ra_mods\", # with ra mods, with feasible conv crit mod\n",
+ "]\n",
+ "\n",
+ "# 3p: triang, psepbqarmax, fimp_14\n",
+ "\n",
+ "qois = [\"w\"]\n",
+ "\n",
+ "# Campaign name bug: string pattern must match unique campaign\n",
+ "campaigns = {}\n",
+ "for campaign_name in campaign_names:\n",
+ " campaigns[campaign_name] = analyse.get_campaign(campaign_name)"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Peek at samples from one campaign: check out a responses."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 105,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# samples = analyse.get_samples(campaigns[\"lt_3p_feas\"])\n",
+ "# samples[qois].describe()"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Distribution of violated constraint residuals"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 106,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " Original convergence criterion | \n",
+ " Feasible convergence criterion | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " 0.057203 | \n",
+ " 0.058435 | \n",
+ "
\n",
+ " \n",
+ " | 1 | \n",
+ " 0.024930 | \n",
+ " 0.062459 | \n",
+ "
\n",
+ " \n",
+ " | 2 | \n",
+ " 0.044947 | \n",
+ " 0.058435 | \n",
+ "
\n",
+ " \n",
+ " | 3 | \n",
+ " 0.084324 | \n",
+ " 0.064453 | \n",
+ "
\n",
+ " \n",
+ " | 4 | \n",
+ " 0.044554 | \n",
+ " 0.068452 | \n",
+ "
\n",
+ " \n",
+ " | ... | \n",
+ " ... | \n",
+ " ... | \n",
+ "
\n",
+ " \n",
+ " | 145 | \n",
+ " 0.038397 | \n",
+ " 0.035472 | \n",
+ "
\n",
+ " \n",
+ " | 146 | \n",
+ " 0.038397 | \n",
+ " 0.034071 | \n",
+ "
\n",
+ " \n",
+ " | 147 | \n",
+ " 0.038397 | \n",
+ " 0.035472 | \n",
+ "
\n",
+ " \n",
+ " | 148 | \n",
+ " 0.025731 | \n",
+ " 0.028744 | \n",
+ "
\n",
+ " \n",
+ " | 149 | \n",
+ " 0.030651 | \n",
+ " 0.028744 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "
150 rows × 2 columns
\n",
+ "