compute inversions for elm

functions/GPT/optimize_parameters.m

-function params = optimize_parameters(measurements, initial_params)
+function params = optimize_parameters(measurements, initial_params, alternative_resolution)
 
-    T = load("data/transfer_matrix.mat");
+    if alternative_resolution
+        T = load("data/transfer_matrix100x40.mat");
+        coords = load("data/coords100x40.mat");
+    else    
+        T = load("data/transfer_matrix.mat");
+        coords = load("data/coords.mat");
+    end
+        
     T = T.transfer_matrix;
-    coords = load("data/coords.mat");
+    
+    for channel = 1:size(T, 1)
+        if measurements(channel) <= 0
+            T(channel, :) = zeros(1, size(T, 2));
+        end
+    end
     
     function value = neg_log_likelihood(params)
         kernel = get_kernel(params, coords.r_coords, coords.z_coords);

functions/GPT/perform_inversion.m

-function post_mean = perform_inversion(measurements, params, transfer_mat, r_coords, z_coords, kernel_inv)
+function post_mean = perform_inversion(measurements, params, transfer_mat, r_coords, z_coords, kernel_inv, alternative_resolution)
 %PERFORM_INVERSION Summary of this function goes here
 %   Detailed explanation goes here
 
@@ -8,6 +8,12 @@ if isempty(kernel_inv)
 end
 variance = (0.0005 * measurements + 0.001).^2;
 
+for channel = 1:size(transfer_mat, 1)
+    if measurements(channel) <= 0
+        transfer_mat(channel, :) = zeros(1, size(transfer_mat, 2));
+    end
+end
+
 
 sigma_d_inv = diag(1./variance);
 
@@ -15,6 +21,10 @@ sigma_post_inv = transfer_mat' * sigma_d_inv * transfer_mat + kernel_inv;
 %cond(sigma_post_inv)
 post_mean = sigma_post_inv \ transfer_mat' * sigma_d_inv * measurements;
 
-post_mean = reshape(post_mean, 25, 75)';
+if alternative_resolution
+    post_mean = reshape(post_mean, 40, 100)';
+else
+    post_mean = reshape(post_mean, 25, 75)';
+end
 end
 

functions/calculate_calibration.m

-function c_factor = calculate_calibration(emissivity, axuv_raw, t_emiss, t_axuv, transfer_matrix)
+function c_factor = calculate_calibration(emissivity, axuv_raw, t_emiss, t_axuv, transfer_matrix, alternative_resolution)
 %CALCULATE_CALIBRATION Summary of this function goes here
 %   Detailed explanation goes here
 
-
 dt_emiss = (t_emiss(101)-t_emiss(1))/100;
 dt_axuv = (t_axuv(101)-t_axuv(1))/100;
 
@@ -12,11 +11,17 @@ axuv_start_index = round((t_emiss(1)-t_axuv(1))/dt_axuv);
 time_multiplier = dt_emiss/dt_axuv;
 c_factor = zeros(size(t_emiss, 1), size(transfer_matrix, 1));
 
+if alternative_resolution
+    num_pixels = 4000;
+else
+    num_pixels = 1875;
+end
+
 for time_index = 1:size(t_emiss, 1)
     emiss = emissivity(:, :, time_index);
     axuv_index = round(time_index*time_multiplier + axuv_start_index);
     axuv = mean(axuv_raw(axuv_index-to_average_over:axuv_index+to_average_over, :));
-    measurement = transfer_matrix * reshape(emiss, 1875, 1);
+    measurement = transfer_matrix * reshape(emiss, num_pixels, 1);
     c_factor(time_index, :) = measurement./axuv';
 end
 

functions/crpptbx-9.2.0/tdi.m

+classdef (HandleCompatible) tdi
+  properties
+    data       = [];
+    dim        = {};
+    units      = '';
+    dimunits   = {};
+    help       = '';
+    validation = [];
+    status     = [];
+    userdata   = [];
+  end
+  
+  methods
+    function y = tdi(arg1,varargin)
+      
+      % TDI class constructor
+      %   A TDI object contains the following fields: _X = EVALUATE(EXPR,P1,...)
+      %          data: DATA(_X)
+      %        dim{i}: DIM_OF(_X,_I-1)
+      %         units: UNITS_OF(_X)
+      %   dimunits{i}: UNITS_OF(DIM_OF(_X,_I-1))
+      %          help: HELP_OF(_X)
+      %    validation: VALIDATION_OF(_X)
+      %        status: return status
+      %      userdata: user use
+      %   TDI(EXPR,P1,...) builds a TDI object using the TDI expression EXPR with Ps
+      %   as $s placeholders. See MDSVALUE for TDI to Matlab and Matlab to TDI
+      %   conversion. X = TDI({FLD1,...},P1,...) fills the FLDs with respective Ps.
+      %   X.FLD refers to FLD. X.FLD = B assigns B to FLD.
+      
+      if nargin > 0 && ~isempty(arg1)
+        
+        % parse arg1 argument
+        switch class(arg1)
+          
+          % tdi(expr,p1,...)
+          case 'char'
+            arg1 = strtrim(arg1);
+            cmd_end = [strfind(arg1,','),strfind(arg1,';')];
+            paren{1} = strfind(arg1,'(');
+            paren{2} = strfind(arg1,')');
+            paren{3} = strfind(arg1,'[');
+            paren{4} = strfind(arg1,']');
+            is = 0;
+            for ii = length(cmd_end):-1:1
+              jj = cmd_end(ii);
+              if jj < length(arg1) && ...
+                  sum(paren{1} < jj) == sum(paren{2} < jj) && ...
+                  sum(paren{3} < jj) == sum(paren{4} < jj)
+                is = jj;
+                break
+              end
+            end
+            % Check if there is any problem with the request
+            %   and reference the result
+            cmd = [arg1(1:is),'__TDIX = ',arg1(is+1:end),',$MISSING'];
+            
+            [y.data,y.status] = mdsdata(cmd,varargin{:});
+            % NOTE: Checks for valid data
+            % 1. Check if status is even
+            % 2. OR status is odd and we have empty dims, and a string in data
+            %       Possible messages are: "No data."
+            %                              "Program XXX not registered"
+            %                              "Program XXX cannot be used for this shot"
+            %               If shot<49830: "Some of the required source nodes to run XXX do not exist"
+            if rem(y.status,2)
+              % Evaluate the result now and store it to avoid re-evaluation
+              % NOTE: This line cannot be mixed with other statements
+              % because of the sequence of COMPILE directives.
+              y.data = mdsdata('__TDIX=`(__TDIX)');
+              pat = {'No data.|',...
+                     'Program (\w+) not registered',...
+                     'Program (\w+) cannot be used for this shot',...
+                     'Some of the required source nodes to run (\w+) do not exist'};
+              if ischar(y.data) && (~isempty(regexp(y.data,pat{1},'once')) || ...
+                                    ~isempty(regexp(y.data,pat{2},'once')) || ...
+                                    ~isempty(regexp(y.data,pat{3},'once')) || ...
+                                    ~isempty(regexp(y.data,pat{4},'once')))
+                y.status=-2;
+              else
+                % NOTE: ndims was processing dimensions from 7 to 0 to find
+                % the first one that existed. ndims2 goes the other way,
+                % from 0 to 7 and looks for the first error. This should be
+                % faster as most signals have 1 or 2 dimensions.
+                % NOTE: If ndims2 returns 0 for the signal, check the data
+                % (if the signal contains a reference to a private
+                % variable, then it will fail to get any dimension).
+                ndim = mdsdata('_n=NDIMS(__TDIX),_n>0 ? _n : NDIMS(DATA(__TDIX))');
+                [y.dim,y.dimunits] = deal(cell(ndim,ndim>0));
+                for k = 1:ndim
+                  y.dim{k}      = mdsdata('__DIM=IF_ERROR(DIM_OF(__TDIX,$1),$MISSING)',k-1);
+                  y.dimunits{k} = mdsdata('IF_ERROR(UNITS_OF(__DIM),"")');
+                end
+                y.units      = mdsdata('IF_ERROR(UNITS_OF(__TDIX),"")');
+                y.help       = mdsdata('IF_ERROR(HELP_OF(__TDIX),"")');
+                y.validation = mdsdata('IF_ERROR(VALIDATION_OF(__TDIX),$ROPRAND)');
+                mdsdata('__TDIX = $MISSING,__DIM=$MISSING');
+              end
+            end
+            
+          % tdi({'fld1',...},p1,...)
+          case 'cell'
+            if ~iscellstr(arg1) || length(arg1) ~= length(varargin)
+              error('Specify compatible field-value lists.')
+            end
+            for k = 1:length(arg1)
+              if isprop(y,arg1{k}), y.(arg1{k}) = varargin{k};
+              else, error('Invalid field name.'), end
+            end
+            
+          % type casting
+          case 'tdi'
+            props = {'data','dim','units','dimunits',...
+              'help','validation','status','userdata'};
+            for ii = 1:numel(props)
+              y.(props{ii}) = arg1.(props{ii});
+            end
+            return
+          otherwise, error('Invalid 1st argument.')
+        end
+        
+      end
+    end
+    
+    disp(x)
+    
+    t = fft(x)
+    
+    plot(x,varargin)
+    
+    varargout = resamp(varargin)
+    
+    n = size(x,ndim)
+    
+    h = surf( x, varargin )
+    
+  end
+end