Primo script Octave nel Script Workspace: fit RC multisorgente

% GX-TXT first Octave Script Workspace scientific workflow.
% Multisource RC discharge alignment, exponential fit, Dataset publication,
% and persistent 2D Plot Dataset graph.
%
% SPDX-License-Identifier: GPL-3.0-or-later

% -------------------------------------------------------------------------
% 1. Inspect the GX-TXT script execution context.
% -------------------------------------------------------------------------
ctx = gxtxt_context();

fprintf('GX-TXT script: %s\n', ctx.script_name);
fprintf('Run ID: %s\n', ctx.run_id);

% This tutorial requires at least two populated Analysis Sources.
if (!isfield(ctx, 'inputs') || !isstruct(ctx.inputs))
  error('RC tutorial: GX-TXT did not provide Analysis Source inputs.');
endif

aliases = fieldnames(ctx.inputs);
source_ids = {};
for k = 1:numel(aliases)
  item = ctx.inputs.(aliases{k});
  if (isstruct(item) && isfield(item, 'dataset_id') && !isempty(item.dataset_id))
    source_ids{end+1} = item.dataset_id;
  endif
endfor

source_ids = unique(source_ids, 'stable');

if (numel(source_ids) < 2)
  error('RC tutorial: select at least two Dataset sources in Analysis.');
endif

fprintf('Input Datasets: %d\n', numel(source_ids));

% -------------------------------------------------------------------------
% 2. Require the public GX-TXT APIs used by the script.
% -------------------------------------------------------------------------
gxtxt_api_require({
  'CANONICAL_DATASET_API_v2@2',
  'DATASET_VARIABLE_SCHEMA_API_v3@3',
  'NUMERIC_API_v1@1:math.numeric',
  'DATASET_PUBLICATION_API_v2@2',
  'PLOT_DATASET_API_v3@3:plot.dataset+plotspec',
  'SESSION_PERSISTENCE_API_v1@1'
});

% -------------------------------------------------------------------------
% 3. Read every source Dataset through public Dataset APIs.
%
% The tutorial expects canonical variables named:
%   time_s
%   voltage_V
% -------------------------------------------------------------------------
time_cells = cell(1, numel(source_ids));
voltage_cells = cell(1, numel(source_ids));
source_variable_ids = {};
time_unit = '';
voltage_unit = '';

for s = 1:numel(source_ids)
  dataset_id = source_ids{s};

  schema = gxtxt_api_call(
      'DATASET_VARIABLE_SCHEMA_API_v3',
      'list',
      {dataset_id});

  time_id = find_variable_id(schema, 'time_s');
  voltage_id = find_variable_id(schema, 'voltage_V');

  if (isempty(time_id) || isempty(voltage_id))
    error('RC tutorial: Dataset %s must contain time_s and voltage_V.', dataset_id);
  endif

  meta = gxtxt_api_call(
      'CANONICAL_DATASET_API_v2',
      'describe',
      {dataset_id});

  block = gxtxt_api_call(
      'CANONICAL_DATASET_API_v2',
      'read_range',
      {dataset_id, {time_id, voltage_id}, 0, meta.row_count});

  t = block.(time_id).values(:);
  v = block.(voltage_id).values(:);

  good = isfinite(t) & isfinite(v);
  t = t(good);
  v = v(good);

  if (numel(t) < 5)
    error('RC tutorial: Dataset %s contains too few finite points.', dataset_id);
  endif

  % NUMERIC_API family alignment requires strictly increasing coordinates.
  [t, order] = sort(t);
  v = v(order);

  if (any(diff(t) <= 0))
    error('RC tutorial: Dataset %s contains duplicate/non-increasing time values.', dataset_id);
  endif

  if (any(v <= 0))
    error('RC tutorial: voltage must remain positive for the log-linear RC fit.');
  endif

  time_cells{s} = t;
  voltage_cells{s} = v;

  tdesc = variable_by_id(schema, time_id);
  vdesc = variable_by_id(schema, voltage_id);

  if (isempty(time_unit)), time_unit = tdesc.unit; endif
  if (isempty(voltage_unit)), voltage_unit = vdesc.unit; endif

  source_variable_ids{end+1} = time_id;
  source_variable_ids{end+1} = voltage_id;

  fprintf('  %s: rows=%d range=[%.6g, %.6g]\n', ...
          dataset_id, numel(t), t(1), t(end));
endfor

source_variable_ids = unique(source_variable_ids, 'stable');

% -------------------------------------------------------------------------
% 4. Build one common time grid.
%
% We use only the overlap shared by every source, so interpolation never
% requires extrapolation.
% -------------------------------------------------------------------------
common_start = -Inf;
common_stop = Inf;

for s = 1:numel(time_cells)
  common_start = max(common_start, time_cells{s}(1));
  common_stop = min(common_stop, time_cells{s}(end));
endfor

if (!(common_stop > common_start))
  error('RC tutorial: the selected sources have no common time interval.');
endif

common_count = 200;
time_common = linspace(common_start, common_stop, common_count).';

% -------------------------------------------------------------------------
% 5. Interpolate all source curves through NUMERIC_API_v1.
%
% The Script Bridge exposes the mathematical-foundation family alignment
% operation.  It accepts typed transport nodes and returns one matrix column
% per source curve.
% -------------------------------------------------------------------------
empty_node = struct('kind', 'empty');
empty_options = struct('kind', 'struct', 'fields', struct());

time_nodes = cell(1, numel(time_cells));
voltage_nodes = cell(1, numel(voltage_cells));

for s = 1:numel(time_cells)
  time_nodes{s} = numeric_array_node(time_cells{s});
  voltage_nodes{s} = numeric_array_node(voltage_cells{s});
endfor

time_family = struct('kind', 'cell', 'items', {time_nodes});
voltage_family = struct('kind', 'cell', 'items', {voltage_nodes});
target_node = numeric_array_node(time_common);

aligned = gxtxt_api_call(
    'NUMERIC_API_v1',
    'foundation',
    {'family', 'align', time_family, voltage_family, target_node, empty_options});

voltage_aligned = aligned.value_real;

if (rows(voltage_aligned) != numel(time_common) || ...
    columns(voltage_aligned) != numel(source_ids))
  error('RC tutorial: unexpected aligned matrix shape.');
endif

voltage_mean = mean(voltage_aligned, 2);

% -------------------------------------------------------------------------
% 6. Fit the zero-offset RC discharge model:
%
%      V(t) = V0 * exp(-t/tau)
%
% Taking the logarithm:
%
%      ln(V) = ln(V0) - t/tau
%
% The logarithm is obtained through NUMERIC_API array_unary, and the
% overdetermined linear system is solved through NUMERIC_API linear_matrix.
% -------------------------------------------------------------------------
positive = isfinite(voltage_mean) & voltage_mean > 0;
t_fit = time_common(positive);
v_fit_input = voltage_mean(positive);

if (numel(t_fit) < 3)
  error('RC tutorial: too few positive consensus points for fitting.');
endif

log_operand = struct('kind', 'vector', 'value', v_fit_input);
log_result = gxtxt_api_call(
    'NUMERIC_API_v1',
    'array_unary',
    {'log', log_operand});
log_voltage = log_result.values(:);

design = [ones(numel(t_fit), 1), t_fit];
rhs = log_voltage;

matrix_A = struct('real', design);
matrix_b = struct('real', rhs);

fit_result = gxtxt_api_call(
    'NUMERIC_API_v1',
    'linear_matrix',
    {'solve', matrix_A, matrix_b});

beta = fit_result.value_real(:);

if (numel(beta) != 2)
  error('RC tutorial: unexpected fit coefficient count.');
endif

log_V0 = beta(1);
slope = beta(2);

if (!(slope < 0))
  error('RC tutorial: fitted slope is not negative; this is not an RC discharge.');
endif

V0_fit = exp(log_V0);
tau_fit = -1 / slope;

fprintf('Fit result: V0 = %.8g V, tau = %.8g s\n', V0_fit, tau_fit);

% -------------------------------------------------------------------------
% 7. Generate the fitted curve through NUMERIC_API array operations.
% -------------------------------------------------------------------------
t_operand = struct('kind', 'vector', 'value', time_common);
slope_operand = struct('kind', 'scalar', 'value', slope);
offset_operand = struct('kind', 'scalar', 'value', log_V0);

scaled_t = gxtxt_api_call(
    'NUMERIC_API_v1',
    'array_binary',
    {'multiply', t_operand, slope_operand});

affine = gxtxt_api_call(
    'NUMERIC_API_v1',
    'array_binary',
    {'add', struct('kind','vector','value',scaled_t.values), offset_operand});

exp_result = gxtxt_api_call(
    'NUMERIC_API_v1',
    'array_unary',
    {'exp', struct('kind','vector','value',affine.values)});

voltage_fit = exp_result.values(:);
residual = voltage_mean - voltage_fit;

rmse = sqrt(mean(residual.^2));

% -------------------------------------------------------------------------
% 8. Publish the derived result as a canonical GX-TXT Dataset.
% -------------------------------------------------------------------------
result_dataset_id = 'tutorial.rc_discharge_fit';

definition = struct();
definition.dataset_id = result_dataset_id;
definition.contract = 'table.generic';
definition.storage = 'auto';
definition.replace_owned = true;
definition.producer_version = '1.0.0';
definition.source_dataset_ids = source_ids;
definition.source_variable_ids = source_variable_ids;
definition.notes = sprintf('Multisource RC discharge fit: tau=%.9g s, V0=%.9g V', ...
                           tau_fit, V0_fit);

writer = gxtxt_api_call(
    'DATASET_PUBLICATION_API_v2',
    'begin',
    {definition});

try
  gxtxt_api_call(
      'DATASET_PUBLICATION_API_v2',
      'declare_dimension',
      {writer, 'row', numel(time_common), struct()});

  time_schema = struct(
      'dtype', 'float64',
      'rank', 1,
      'dimension_ids', 'row',
      'coordinate_role', 'coordinate',
      'semantic_role', 'independent',
      'unit', time_unit);

  mean_schema = struct(
      'dtype', 'float64',
      'rank', 1,
      'dimension_ids', 'row',
      'coordinate_role', 'data',
      'semantic_role', 'measured',
      'unit', voltage_unit);

  fit_schema = struct(
      'dtype', 'float64',
      'rank', 1,
      'dimension_ids', 'row',
      'coordinate_role', 'data',
      'semantic_role', 'derived',
      'unit', voltage_unit);

  residual_schema = fit_schema;

  gxtxt_api_call(
      'DATASET_PUBLICATION_API_v2',
      'declare_variable',
      {writer, 'v_time', time_schema});

  gxtxt_api_call(
      'DATASET_PUBLICATION_API_v2',
      'declare_variable',
      {writer, 'v_mean', mean_schema});

  gxtxt_api_call(
      'DATASET_PUBLICATION_API_v2',
      'declare_variable',
      {writer, 'v_fit', fit_schema});

  gxtxt_api_call(
      'DATASET_PUBLICATION_API_v2',
      'declare_variable',
      {writer, 'v_residual', residual_schema});

  gxtxt_api_call(
      'DATASET_PUBLICATION_API_v2',
      'write_variable',
      {writer, 'v_time', time_common, numel(time_common)});

  gxtxt_api_call(
      'DATASET_PUBLICATION_API_v2',
      'write_variable',
      {writer, 'v_mean', voltage_mean, numel(voltage_mean)});

  gxtxt_api_call(
      'DATASET_PUBLICATION_API_v2',
      'write_variable',
      {writer, 'v_fit', voltage_fit, numel(voltage_fit)});

  gxtxt_api_call(
      'DATASET_PUBLICATION_API_v2',
      'write_variable',
      {writer, 'v_residual', residual, numel(residual)});

  published = gxtxt_api_call(
      'DATASET_PUBLICATION_API_v2',
      'finalize',
      {writer});

  writer = '';
catch err
  if (!isempty(writer))
    gxtxt_api_call(
        'DATASET_PUBLICATION_API_v2',
        'abort',
        {writer});
  endif
  rethrow(err);
end

fprintf('Published Dataset: %s\n', published.dataset_id);

% -------------------------------------------------------------------------
% 9. Create one simple persistent GX-TXT 2D graph.
% -------------------------------------------------------------------------
plot_spec = gxtxt_api_call(
    'PLOT_DATASET_API_v3',
    'new_spec',
    {'plot.dataset.xy', 'RC discharge: multisource mean and exponential fit'});

panel = plot_spec.panels{1};

measured_series = struct(
    'id', 'mean_measured',
    'dataset_id', result_dataset_id,
    'x_variable_id', 'v_time',
    'y_variable_id', 'v_mean',
    'axis', 'left',
    'kind', 'line_marker',
    'label', 'Mean measured');

fit_series = struct(
    'id', 'exponential_fit',
    'dataset_id', result_dataset_id,
    'x_variable_id', 'v_time',
    'y_variable_id', 'v_fit',
    'axis', 'left',
    'kind', 'line',
    'label', sprintf('Exponential fit, tau = %.4g s', tau_fit));

panel.series = {measured_series, fit_series};
plot_spec.panels = {panel};

gxtxt_api_call(
    'PLOT_DATASET_API_v3',
    'validate_spec',
    {plot_spec});

plot_ref = gxtxt_api_call(
    'PLOT_DATASET_API_v3',
    'upsert_plot',
    {'tutorial_rc_discharge_fit', plot_spec});

gxtxt_api_call(
    'PLOT_DATASET_API_v3',
    'render',
    {plot_ref});

% -------------------------------------------------------------------------
% 10. Save a small human-readable report in this script run directory.
% -------------------------------------------------------------------------
report_path = fullfile(ctx.output_dir, 'rc_fit_report.txt');
fid = fopen(report_path, 'w');
if (fid >= 0)
  fprintf(fid, 'GX-TXT RC discharge tutorial\n');
  fprintf(fid, '============================\n\n');
  fprintf(fid, 'Input Dataset count: %d\n', numel(source_ids));
  fprintf(fid, 'Common interval: %.9g .. %.9g %s\n', ...
          common_start, common_stop, time_unit);
  fprintf(fid, 'Common grid points: %d\n', numel(time_common));
  fprintf(fid, 'Fitted V0: %.12g %s\n', V0_fit, voltage_unit);
  fprintf(fid, 'Fitted tau: %.12g %s\n', tau_fit, time_unit);
  fprintf(fid, 'Consensus RMSE: %.12g %s\n', rmse, voltage_unit);
  fprintf(fid, 'Published Dataset: %s\n', result_dataset_id);
  fclose(fid);
endif

% Persist the newly published Dataset and persistent Plot Dataset Variant.
gxtxt_api_call(
    'SESSION_PERSISTENCE_API_v1',
    'commit',
    {});

fprintf('PASS RC tutorial: Dataset=%s tau=%.8g %s variant=%s\n', ...
        result_dataset_id, tau_fit, time_unit, plot_ref.variant_id);

% =========================================================================
% Local helper functions
% =========================================================================

function variable_id = find_variable_id(schema, canonical_name)
  variable_id = '';
  for i = 1:numel(schema)
    item = schema{i};
    if (isstruct(item) && isfield(item, 'canonical') && ...
        strcmp(item.canonical, canonical_name))
      variable_id = item.variable_id;
      return;
    endif
  endfor
endfunction

function item = variable_by_id(schema, variable_id)
  item = struct();
  for i = 1:numel(schema)
    candidate = schema{i};
    if (isstruct(candidate) && isfield(candidate, 'variable_id') && ...
        strcmp(candidate.variable_id, variable_id))
      item = candidate;
      return;
    endif
  endfor
  error('RC tutorial: Variable ID not found in schema: %s', variable_id);
endfunction

function node = numeric_array_node(values)
  values = values(:);
  node = struct(
      'kind', 'array',
      'shape', [numel(values), 1],
      'real', values.');
endfunction