text:
function [contour] = dicomrt_plotVOI(dir,slice,VOI,voi2use,xmesh,ymesh,zmesh,PatientPosition)
% dicomrt_plorVOI(dir,slice,VOI,voi2use,xmesh,ymesh,zmesh,PatientPosition)
%
% Plot VOIs in X and Y plane.
%
% dir is the direction (=1 -> x, =2 -> y, =3 -> z)
% slice is the slice number to plot the contour for
% VOI is a cell array which contain the patients VOIs.
% voi2use is a vector pointing to the number of VOIs to be displayed
%
% Differs from dicomrt_overlayvoi as voi section to plot is not searched.
%
% See also dicomrt_findsliceVECT, dicomrt_displaycontour, dicomrt_loadvoi, dicomrt_rendervoi
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Check number of argument and set-up some parameters and variables
error(nargchk(8,9,nargin))
[VOI_temp]=dicomrt_checkinput(VOI);
VOI=dicomrt_varfilter(VOI_temp);

voitype=dicomrt_checkvoitype(VOI_temp);

% Set parameter
voilookup=1;
voiZ=[];
C=[];
h=[];

% Initialize plot parameters:
color = char('r','g','b','c','m','y','k','w');
line = char('-','--',':','-.');
width = [0.5,1.0,1.5,2.0];

% Default values
colorselect=8;
lineselect=1;
linewidth=2;

if dir==1 & isequal(voitype,'3D')==1
    for i=1:length(voi2use)
        dicomrt_plotcontourc(VOI{voi2use(i),2}{1}{slice}{2});
    end
elseif dir==1 & isequal(voitype,'3D')~=1
    warning('Contour plot is not possible. Use dicomrt_build3dVOI to build a 3D VOI');
elseif dir==2 & isequal(voitype,'3D')==1
    for i=1:length(voi2use)
        dicomrt_plotcontourc(VOI{voi2use(i),2}{2}{slice}{2});
    end
elseif dir==2 & isequal(voitype,'3D')~=1
    warning('Contour plot is not possible. Use dicomrt_build3dVOI to build a 3D VOI');
elseif dir==3 & isequal(voitype,'3D')~=1
    for i=1:length(voi2use)
        [voiZ,voiZ_index]=dicomrt_getvoiz(VOI_temp,voi2use(i));
        voiZ=dicomrt_makevertical(voiZ);
        zslice=dicomrt_findpointVECT(voiZ,zmesh(slice));
        zslice=voiZ_index(zslice);
        if isempty(zslice)~=1
            closexVOI=[VOI{voi2use(i),2}{zslice}(:,1);VOI{voi2use(i),2}{zslice}(1,1)];
            closeyVOI=[VOI{voi2use(i),2}{zslice}(:,2);VOI{voi2use(i),2}{zslice}(1,2)];
            hplot=plot(closexVOI,closeyVOI);
            set(hplot,'Color',color(colorselect,:));
            set(hplot,'LineStyle',line(lineselect,:));
            set(hplot,'LineWidth', linewidth);
        else
            warning(['dicomrt_plot3dVOI: VOI :',VOI{voi2use(i),1},' is not defined in the selected slice']);
        end
    end
elseif dir==3 & isequal(voitype,'3D')==1
    for i=1:length(voi2use)
        [voiZ,voiZ_index]=dicomrt_getvoiz(VOI_temp,voi2use(i));
        voiZ=dicomrt_makevertical(voiZ);
        zslice=dicomrt_findpointVECT(voiZ,zmesh(slice));
        zslice=voiZ_index(zslice);
        if isempty(zslice)~=1
            closexVOI=[VOI{voi2use(i),2}{3}{zslice}(:,1);VOI{voi2use(i),2}{3}{zslice}(1,1)];
            closeyVOI=[VOI{voi2use(i),2}{3}{zslice}(:,2);VOI{voi2use(i),2}{3}{zslice}(1,2)];
            hplot=plot(closexVOI,closeyVOI);
            set(hplot,'Color',color(colorselect,:));
            set(hplot,'LineStyle',line(lineselect,:));
            set(hplot,'LineWidth', linewidth);
        else
            warning(['dicomrt_plot3dVOI: VOI :',VOI{voi2use(i),1},' is not defined in the selected slice']);
        end
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ a, b ] = p26_lim ( )

%*****************************************************************************80
%
%% P26_LIM returns the integration limits for problem 26.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 2.0 * pi;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ r, lchar, ierror ] = s_to_r8 ( s )

%*****************************************************************************80
%
%% S_TO_R8 reads an R8 from a string.
%
%  Discussion:
%
%    This routine will read as many characters as possible until it reaches
%    the end of the string, or encounters a character which cannot be
%    part of the real number.
%
%    Legal input is:
%
%       1 blanks,
%       2 '+' or '-' sign,
%       2.5 spaces
%       3 integer part,
%       4 decimal point,
%       5 fraction part,
%       6 'E' or 'e' or 'D' or 'd', exponent marker,
%       7 exponent sign,
%       8 exponent integer part,
%       9 exponent decimal point,
%      10 exponent fraction part,
%      11 blanks,
%      12 final comma or semicolon.
%
%    with most quantities optional.
%
%  Example:
%
%    S                 R
%
%    '1'               1.0
%    '     1   '       1.0
%    '1A'              1.0
%    '12,34,56'        12.0
%    '  34 7'          34.0
%    '-1E2ABCD'        -100.0
%    '-1X2ABCD'        -1.0
%    ' 2E-1'           0.2
%    '23.45'           23.45
%    '-4.2E+2'         -420.0
%    '17d2'            1700.0
%    '-14e-2'         -0.14
%    'e2'              100.0
%    '-12.73e-9.23'   -12.73 * 10.0**(-9.23)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string containing the
%    data to be read.  Reading will begin at position 1 and
%    terminate at the end of the string, or when no more
%    characters can be read to form a legal real.  Blanks,
%    commas, or other nonnumeric data will, in particular,
%    cause the conversion to halt.
%
%    Output, real R, the value that was read from the string.
%
%    Output, integer LCHAR, the number of characters of S that were used to form R.
%
%    Output, integer IERROR, is 0 if no error occurred.
%
  s_length = s_len_trim ( s );
  ierror = 0;
  r = 0.0;
  lchar = -1;
  isgn = 1;
  rtop = 0.0;
  rbot = 1.0;
  jsgn = 1;
  jtop = 0;
  jbot = 1;
  ihave = 1;
  iterm = 0;

  while ( 1 )

    lchar = lchar + 1;
    c = s(lchar+1);
%
%  Blank character.
%
    if ( c == ' ' )

      if ( ihave == 2 )

      elseif ( ihave == 6 | ihave == 7 )
        iterm = 1;
      elseif ( 1 < ihave )
        ihave = 11;
      end
%
%  Comma.
%
    elseif ( c == ',' | c == ';' )

      if ( ihave ~= 1 )
        iterm = 1;
        ihave = 12;
        lchar = lchar + 1;
      end
%
%  Minus sign.
%
    elseif ( c == '-' )

      if ( ihave == 1 );
        ihave = 2;
        isgn = -1;
      elseif ( ihave == 6 )
        ihave = 7;
        jsgn = -1;
      else
        iterm = 1;
      end
%
%  Plus sign.
%
    elseif ( c == '+' )

      if ( ihave == 1 )
        ihave = 2;
      elseif ( ihave == 6 )
        ihave = 7;
      else
        iterm = 1;
      end
%
%  Decimal point.
%
    elseif ( c == '.' )

      if ( ihave < 4 )
        ihave = 4;
      elseif ( 6 <= ihave & ihave <= 8 )
        ihave = 9;
      else
        iterm = 1;
      end
%
%  Exponent marker.
%
    elseif ( ch_eqi ( c, 'E' ) | ch_eqi ( c, 'D' ) )

      if ( ihave < 6 )
        ihave = 6;
      else
        iterm = 1;
      end
%
%  Digit.
%
    elseif ( ihave < 11 & ch_is_digit ( c ) )

      if ( ihave <= 2 )
        ihave = 3;
      elseif ( ihave == 4 )
        ihave = 5;
      elseif ( ihave == 6 | ihave == 7 )
        ihave = 8;
      elseif ( ihave == 9 )
        ihave = 10;
      end

      d = ch_to_digit ( c );

      if ( ihave == 3 )
        rtop = 10.0 * rtop + d;
      elseif ( ihave == 5 )
        rtop = 10.0 * rtop + d;
        rbot = 10.0 * rbot;
      elseif ( ihave == 8 )
        jtop = 10 * jtop + d;
      elseif ( ihave == 10 )
        jtop = 10 * jtop + d;
        jbot = 10 * jbot;
      end
%
%  Anything else is regarded as a terminator.
%
    else
      iterm = 1;
    end
%
%  If we haven't seen a terminator, and we haven't examined the
%  entire string, go get the next character.
%
    if ( iterm == 1 | s_length <= lchar + 1 )
      break;
    end

  end
%
%  If we haven't seen a terminator, and we have examined the
%  entire string, then we're done, and LCHAR is equal to S_LENGTH.
%
  if ( iterm ~= 1 & lchar + 1 == s_length )
    lchar = s_length;
  end
%
%  Number seems to have terminated.  Have we got a legal number?
%  Not if we terminated in states 1, 2, 6 or 7!
%
  if ( ihave == 1 | ihave == 2 | ihave == 6 | ihave == 7 )
    ierror = ihave;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'S_TO_R8 - Fatal error!\n' );
    fprintf ( 1, '  IHAVE = %d\n', ihave );
    error ( 'S_TO_R8 - Fatal error!' );
  end
%
%  Number seems OK.  Form it.
%
  if ( jtop == 0 )
    rexp = 1.0;
  else

    if ( jbot == 1 )
      rexp = 10.0^( jsgn * jtop );
    else
      rexp = jsgn * jtop;
      rexp = rexp / jbot;
      rexp = 10.0^rexp;
    end

  end

  r = isgn * rexp * rtop / rbot;

  return
end


text:
function res = spm_eeg_reduce_imagcsd(S)
% Plugin for data reduction based on the imaginary part of CSD
% with a reference chhannel
% FORMAT res = spm_eeg_reduce_imagcsd(S)
%
% S                     - input structure
% fields of S:
%
%
% Output:
%  res -
%   If no input is provided the plugin returns a cfg branch for itself
%
%   If input is provided:
%      montage struct implementing projection to PCA subspace
%______________________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_reduce_imagcsd.m 6852 2016-08-01 12:48:44Z vladimir $

if nargin == 0

    origchan = cfg_branch;
    origchan.tag = 'origchan';
    origchan.name = 'Channels to reduce';
    origchan.val = {spm_cfg_eeg_channel_selector};

    refchan = cfg_branch;
    refchan.tag = 'refchan';
    refchan.name = 'Reference channels';
    refchan.val = {spm_cfg_eeg_channel_selector};

    outlabel = cfg_entry;
    outlabel.tag = 'outlabel';
    outlabel.name = 'Output channel label';
    outlabel.strtype = 's';
    outlabel.num = [1 Inf];
    outlabel.help = {'Label for the output channel(s).'};

    foi = cfg_entry;
    foi.tag = 'foi';
    foi.name = 'Frequency band of interest';
    foi.strtype = 'r';
    foi.num = [1 2];
    foi.val = {[0 Inf]};
    foi.help = {'Frequency window to optimize for'};

    chanset = cfg_branch;
    chanset.tag = 'chanset';
    chanset.name = 'Set';
    chanset.val = {origchan, refchan, outlabel, foi};

    chansets = cfg_repeat;
    chansets.tag = 'chansets';
    chansets.name = 'Channel sets';
    chansets.values = {chanset};
    chansets.num = [1 Inf];
    chansets.val = {chanset};

    imagcsd = cfg_branch;
    imagcsd.tag = 'imagcsd';
    imagcsd.name = 'Imaginary part of CSD';
    imagcsd.val = {chansets};

    res = imagcsd;

    return
end

flag_tbx = license('checkout','signal_toolbox') && ~isempty(ver('signal'));
if flag_tbx
    taper = 'dpss';
else
    taper = 'sine';
end

D = S.D;

nsets = numel(S.chanset);
badind = D.badchannels;

% Assuming projecting to columns
montage = [];
montage.labelorg    = D.chanlabels;
montage.labelnew    = {};
montage.chantypenew = {};
montage.tra         = zeros(0, D.nchannels);

spm('Pointer', 'Watch');drawnow;
spm_progress_bar('Init', nsets, 'Channel sets processed'); drawnow;
if nsets > 100, Ibar = floor(linspace(1, nsets,100));
else Ibar = 1:nsets; end

for i = 1:nsets

    spm_progress_bar('Set','ylabel','preparing data...');

    origind = setdiff(D.selectchannels(spm_cfg_eeg_channel_selector(S.chanset(i).origchan.channels)), badind);
    refind = setdiff(D.selectchannels(spm_cfg_eeg_channel_selector(S.chanset(i).refchan.channels)), badind);

    if length(refind)~=1
        error(sprintf('One reference chhannel is necessary, found %d.', length(refind)));
    end

    Y  = D(origind, :, :);
    Yr = spm_squeeze(D(refind, :, :), 1)';

    if isequal(D.type, 'continuous')
        fs   = floor(D.fsample);
        ind = repmat((1:fs:D.nsamples)', 1, fs);
        ind = ind + repmat(0:(fs-1), size(ind, 1), 1);

        Yr = Yr(ind);

        time = (0:(fs-1))/fs;
    else
        time = D.time;
    end

    spm_progress_bar('Set','ylabel','computing CSD...');

    foi =  S.chanset(i).foi;
    foi(isinf(foi)) = 0.5*D.fsample;

    fY =[];
    for c = 1:size(Y, 1)
        Yc = spm_squeeze(Y(c, :, :), 1)';
        if size(Yc, 1) == 1
            Yc = Yc(ind);
        end
        [spectrum,ntaper,freqoi] = ft_specest_mtmfft(Yc, time, 'taper', taper, 'freqoi', mean(foi), 'tapsmofrq', 0.5*diff(foi), 'verbose', 0);
        fY = [fY spectrum(:)];
    end

    fYr = ft_specest_mtmfft(Yr, time, 'taper', taper, 'freqoi', mean(foi), 'tapsmofrq', 0.5*diff(foi), 'verbose', 0);

    fYr = fYr(:);

    csd = mean(fY.*repmat(conj(fYr), 1, size(fY, 2)));

    mom = imag(csd)./norm(imag(csd));

    montage.labelnew{end+1, 1} = S.chanset(i).outlabel;

    montage.tra(end+1, end)  = 0;
    montage.tra(end, origind) = mom;

    montage.chantypenew{end+1}='LFP';

    if ismember(i, Ibar)
        spm_progress_bar('Set', i); drawnow;
    end
end

spm_progress_bar('Clear');

if ~isempty(S.chanind)
    montage.labelnew = [montage.labelnew; D.chanlabels(S.chanind)'];
    I = eye(D.nchannels);
    montage.tra = [montage.tra; I(S.chanind, :)];
    montage.chantypenew = [montage.chantypenew, D.chantype(S.chanind)];
end

res = montage;

author: spm


text:
function [xn] = normalize_pixel_fisheye(x_kk,fc,cc,kc,alpha_c)

%normalize
%
%[xn] = normalize_pixel(x_kk,fc,cc,kc,alpha_c)
%
%Computes the normalized coordinates xn given the pixel coordinates x_kk
%and the intrinsic camera parameters fc, cc and kc.
%
%INPUT: x_kk: Feature locations on the images
%       fc: Camera focal length
%       cc: Principal point coordinates
%       kc: Fisheye distortion coefficients
%       alpha_c: Skew coefficient
%
%OUTPUT: xn: Normalized feature locations on the image plane (a 2XN matrix)

if nargin < 5,
   alpha_c = 0;
   if nargin < 4;
      kc = [0;0;0;0;0];
      if nargin < 3;
         cc = [0;0];
         if nargin < 2,
            fc = [1;1];
         end;
      end;
   end;
end;

% First: Subtract principal point, and divide by the focal length:
x_distort = [(x_kk(1,:) - cc(1))/fc(1);(x_kk(2,:) - cc(2))/fc(2)];

% Second: undo skew
x_distort(1,:) = x_distort(1,:) - alpha_c * x_distort(2,:);

% Third: Compensate for lens distortion:
xn = comp_fisheye_distortion(x_distort,kc);

repo: ekfmonoslam


text:
function f = neglogsigmoid_fh(w)
% This is a factory for a function handle to an MV2DF, which represents
% the vectorization of the logsigmoid function. The mapping is, in
% MATLAB-style code:
%
%   y = log(sigmoid(w)) = log(1./1+exp(-w)) = -log(1+exp(-w))
%
% Inputs:
%   m: the number of inputs to each individual logsumexp calculation.
%   direction: 1 sums down columns, or 2 sums accross rows.
%   w: optional, if ssupplied
%
% Outputs:
%   f: a function handle to the MV2DF described above.
%
% see: MV2DF_API_DEFINITION.readme

if nargin==0
    test_this();
    return;
end

f = vectorized_function([],@(x)F0(x));

if exist('w','var') && ~isempty(w)
    f = f(w);
end

end

function [y,f1] = F0(x)
logp1 = -neglogsigmoid(x);
logp2 = -neglogsigmoid(-x);
y = -logp1;
f1 = @() F1(logp1,logp2);
end

function [J,f2,linear] = F1(logp1,logp2)
linear = false;
J = -exp(logp2);
f2 = @(dx) F2(dx,logp1,logp2);
end

function h = F2(dx,logp1,logp2)
h = dx.*exp(logp1+logp2);
end

function test_this()
n = 10;
f = neglogsigmoid_fh([]);
x = randn(n,1);
test_MV2DF(f,x);
end

author: nesl


text:
function [E_map,w,f,costf,iter,value] = pixelid(data,F_fw,te,varargin)

% in
%   data[n,n,nsets]     source data
%   F_fw = fat water shift (omega)
%   te [nsets]   time
%options
%   niter  # of iterations
%   x       [array or one value]
%   y       [array or one value]
%           x,y are array or one value of object field map
%           If x = 100, y = 200 then we choose just (100,200)th pixel of field map
%           If x = 1:100, y = 200:202 then we choose (1:100,200:202) of
%           field map ( = estimate 100*103 pixels)
%   init  [n,n] Initial guess of field map
%   cost    If cost = 0 : Do not calculate cost function (Default)
%           If cost = 1 : Calculate cost function
%% initialize

n = size(data,1); % length of the image

%defaults
arg.niter = 30;
arg.x = 1:512;
arg.y = 1:512;  % Just choose (100,100)th pixel of field map
arg.init = zeros(n,n);
arg.cost = 0;
arg = vararg_pair(arg, varargin);

c = cos(F_fw*te).';
d = sin(F_fw*te).'; % for Appendix A and B
A1 =  [1 0 c(1) -d(1); 1 0 c(2) -d(2); 1 0 c(3) -d(3); 0 1 d(1) c(1); 0 1 d(2) c(2); 0 1 d(3) c(3)];
A2 = inv(A1'*A1)*A1';  % for Appendix A
map = arg.init; % initial field map
rgmap = arg.init;
A3(:,1) = ones(length(te),1);
A3(:,2) = exp(j*F_fw.*te); % for calculate cost function
if arg.cost == 1
    length_cost = size(arg.x,2)*size(arg.y,2);
    costf = zeros(length_cost*2,arg.niter*2);
end
iteration = 0;

%% pixel independent method

for arr = arg.x
    if (rem(arr,100) == 0)
         disp(['# of iteration of RE = ' num2str(arr)])
    end

    for col = arg.y

        iteration = iteration+1;
        init_guess = arg.init(arr,col);
        [est_map,cost,iter] = regionpixel(data,arr,col,init_guess,te,c,d,A1,A2,A3,arg.niter,arg.cost);
        if arg.cost == 1
                costf((iteration*2-1):iteration*2,:) = cost;
                % cost(2*k-1,:) = field map values(record of estimated field map values in each iteration)
                % cost(2*k,:) = cost values
        end

        % My constraint
        if est_map > 90
            est_map = 90;
        elseif est_map <-80
            est_map = -80;
        end
        map((col-1)*n+arr) = est_map;
    end
end
if arg.cost == 1
    %% plot the cost functions (plot only the first pixel's cost function)
    % just plot costfunction(field map,water^,fat^),water^,fat^ is estimated
    % water, fat based on field map
    figure(5);
    plot(costf(1,(1:arg.niter)*2-1),costf(2,(1:arg.niter)*2-1),'o');title('cost function');
    hold on
    plot(costf(1,iter*2-1),costf(2,iter*2-1),'o','color','red'); % final estimated field map
    hold off

    %% plot the whole cost function (from -150Hz - 150Hz)
    a = linspace(-1200,1200,10000);
    orig_s = reshape(data(arg.x(1),arg.y(1),:),3,1);
    for k = 1:10000
        s =  orig_s.*(exp(-j*2*pi*a(k)*te).'); % shat (m is freq)
        shat = [real(s);imag(s)];
        rowhat = A2*shat;
        water = rowhat(1) + j*rowhat(2);
        fat = rowhat(3) + j*rowhat(4);
        temp(1,k) = norm(orig_s - A3*[water;fat].*(exp(j*2*pi*a(k)*te).'));
    end
    figure(6);
    plot(a,temp,'o');
    w = 0;
    f = 0;
    E_map = map;
    value = est_map;
else
    costf = 0;
    value = 0;
    E_map = zeros(n,n);
    E_map(find(map>0)) = map(find(map>0));
    E_map(find(map<0)) = map(find(map<0));
    E_map = lpass(E_map);
    E_map = E_map*2*pi;
    rgmap = rgmap*2*pi;
   [w f] = restore(data,E_map,F_fw,te);
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
classdef casadi_block < matlab.System & matlab.system.mixin.Propagates
    % Casadi Implementation of the Momentum Based Velocity COntrol.

    properties
        % Public, tunable properties.

    end

    properties (DiscreteState)
    end

    properties (Access = private)
        casadi_optimizer;
        M;
        h;
        Jc;
        Jcmm;
        nu_k;
        C_friction;
        b_friction;
        s_dot_star;
        H_star;
        tau_k;
        s_dot_k_1;
        v_b_k_1;
        f_k;
        Jc_dot_nu;
        sol;
        tau_meas;
        f_meas;
        solver_fails_counter;
        s;
        CoM_measured;
        K_friction;
        t_step;
        epsilon_CoM;
        CoM_max;
        CoM_min;
        CoM_vel_max;
        CoM_vel_min;
        ActivateComHeightConstraint;
        CoM_Z_reference_velocity;
        feet_correction;
        holonomicConstraintActivation;
        JTorso;
        SaturationMinJoints;
        SaturationMaxJoints;
        J_com;

    end

    methods (Access = protected)
        function num = getNumInputsImpl(~)
            num = 36;
        end
        function num = getNumOutputsImpl(~)
            num = 4;
        end
        function [dt1,dt2,dt3,dt4] = getOutputDataTypeImpl(~)
            dt1 = 'double';
            dt2 = 'double';
            dt3 = 'double';
            dt4 = 'double';

        end

        function [dt1,dt2,dt3,dt4,dt5,dt6,dt7,dt8,dt9,dt10,dt11,dt12,dt13,dt14,dt15,dt16,dt17,dt18,dt19,dt20,dt21,dt22,dt23,dt24, ...
                dt25,dt26,dt27,dt28,dt29,dt30,dt31,dt32,dt33,dt34,dt35,dt36] = getInputDataTypeImpl(~)
            dt1 = 'double';
            dt2 = 'double';
            dt3 = 'double';
            dt4 = 'double';
            dt5 = 'double';
            dt6 = 'double';
            dt7 = 'double';
            dt8 = 'double';
            dt9 = 'double';
            dt10 = 'double';
            dt11 = 'double';
            dt12 = 'double';
            dt13 = 'double';
            dt14 = 'double';
            dt15 = 'double';
            dt16 = 'double';
            dt17 = 'double';
            dt18 = 'double';
            dt19 = 'double';
            dt20 = 'double';
            dt21 = 'double';
            dt22 = 'double';
            dt23 = 'double';
            dt24 = 'double';
            dt25 = 'double';
            dt26 = 'double';
            dt27 = 'double';
            dt28 = 'double';
            dt29 = 'double';
            dt30 = 'double';
            dt31 = 'double';
            dt32 = 'double';
            dt33 = 'double';
            dt34 = 'double';
            dt35 = 'double';
            dt36 = 'double';
        end
        function [sz1,sz2,sz3,sz4] = getOutputSizeImpl(obj)
            % Hack for avoiding hard coding the DoF of the Robot
            size_temp = propagatedInputSize(obj,1);
            sz1 = [6,1];
            sz2 = [size_temp(1)-6,1];
            sz3 = [size_temp(1)-6,1];
            sz4 = [12,1];
        end
        function [sz1,sz2,sz3,sz4,sz5,sz6,sz7,sz8,sz9,sz10,sz11,sz12,sz13,sz14,sz15,sz16, sz17,sz18,sz19,sz20,sz21,sz22, ...
                sz23,sz24,sz25,sz26,sz27,sz28,sz29,sz30,sz31,sz32,sz33,sz34, sz35,sz36] = getInputSizeImpl(~)
            sz1 = [1,1];
            sz2 = [1,1];
            sz3 = [1,1];
            sz4 = [1,1];
            sz5 = [1,1];
            sz6 = [1,1];
            sz7 = [1,1];
            sz8 = [1,1];
            sz9 = [1,1];
            sz10 = [1,1];
            sz11 = [1,1];
            sz12 = [1,1];
            sz13 = [1,1];
            sz14 = [1,1];
            sz15 = [1,1];
            sz16 = [1,1];
            sz17 = [1,1];
            sz18 = [1,1];
            sz19 = [1,1];
            sz20 = [1,1];
            sz21 = [1,1];
            sz22 = [1,1];
            sz23 = [1,1];
            sz24 = [1,1];
            sz25 = [1,1];
            sz26 = [1,1];
            sz27 = [1,1];
            sz28 = [1,1];
            sz29 = [1,1];
            sz30 = [1,1];
            sz31 = [1,1];
            sz32 = [1,1];
            sz33 = [1,1];
            sz34 = [1,1];
            sz35 = [1,1];
            sz36 = [1,1];
        end
        function [cp1,cp2,cp3,cp4,cp5,cp6,cp7,cp8,cp9,cp10,cp11,cp12,cp13, cp14, cp15,cp16, cp17,cp18,cp19,cp20,cp21,...
                cp22,cp23,cp24,cp25,cp26,cp27,cp28,cp29,cp30,cp31,cp32,cp33,cp34,cp35,cp36] = isInputComplexImpl(~)
            cp1 = false;
            cp2 = false;
            cp3 = false;
            cp4 = false;
            cp5 = false;
            cp6 = false;
            cp7 = false;
            cp8 = false;
            cp9 = false;
            cp10 = false;
            cp11 = false;
            cp12 = false;
            cp13 = false;
            cp14 = false;
            cp15 = false;
            cp16 = false;
            cp17 = false;
            cp18 = false;
            cp19 = false;
            cp20 = false;
            cp21 = false;
            cp22 = false;
            cp23 = false;
            cp24 = false;
            cp25 = false;
            cp26 = false;
            cp27 = false;
            cp28 = false;
            cp29 = false;
            cp30 = false;
            cp31 = false;
            cp32 = false;
            cp33 = false;
            cp34 = false;
            cp35 = false;
            cp36 = false;
        end
        function [cp1,cp2,cp3,cp4] = isOutputComplexImpl(~)
            cp1 = false;
            cp2 = false;
            cp3 = false;
            cp4 = false;
        end
        function [fz1,fz2,fz3,fz4,fz5,fz6,fz7,fz8,fz9,fz10, fz11,fz12,fz13,fz14,fz15,fz16, fz17,fz18,fz19,fz20,fz21, ...
                fz22,fz23,fz24,fz25,fz26,fz27,fz28,fz29,fz30,fz31,fz32,fz33,fz34, fz35,fz36] = isInputFixedSizeImpl(~)
            fz1 = true;
            fz2 = true;
            fz3 = true;
            fz4 = true;
            fz5 = true;
            fz6 = true;
            fz7 = true;
            fz8 = true;
            fz9 = true;
            fz10 = true;
            fz11 = true;
            fz12 = true;
            fz13 = true;
            fz14 = true;
            fz15 = true;
            fz16 = true;
            fz17 = true;
            fz18 = true;
            fz19 = true;
            fz20 = true;
            fz21 = true;
            fz22 = true;
            fz23 = true;
            fz24 = true;
            fz25 = true;
            fz26 = true;
            fz27 = true;
            fz28 = true;
            fz29 = true;
            fz30 = true;
            fz31 = true;
            fz32 = true;
            fz33 = true;
            fz34 = true;
            fz35 = true;
            fz36 = true;
        end
        function [fz1,fz2,fz3,fz4] = isOutputFixedSizeImpl(~)
            fz1 = true;
            fz2 = true;
            fz3 = true;
            fz4 = true;
        end
        function setupImpl(obj,M,h,Jc,Jcmm,nu_k,Jc_dot_nu,C_friction,b_friction,H_star,s_dot_star,tau_meas, contact_status, f_meas, x_com, J_com, w_H_b,...
                s,CoM_measured,ADD_FRICTION, K_viscous_friction, T, Gamma, baseVelocitySat, jointsVelocitySat, torqueSat, wrenchesSat, tStep, NDOF, epsilon_CoM, CoM_limits,CoM_vel_limits, activateCoMHeightTask, CoMZReferenceVelocity, correctionFeet, ActivateHolonomicConstraint, JTorso)

            NDOF = size(K_viscous_friction,1);
            obj.solver_fails_counter = 0;
            import casadi.*
            obj.solver_fails_counter = 0;
            optimization_type = 'conic';
            solver = 'osqp';
            p_opts = struct('expand',true);
            s_opts = struct();

            obj.casadi_optimizer = casadi.Opti(optimization_type);
            obj.casadi_optimizer.solver(solver, p_opts, s_opts);

            %Definition of the decision variables
            obj.tau_k = obj.casadi_optimizer.variable(NDOF);
            obj.s_dot_k_1 = obj.casadi_optimizer.variable(NDOF);
            obj.v_b_k_1 = obj.casadi_optimizer.variable(6);
            obj.f_k = obj.casadi_optimizer.variable(12);
            nu_k_1 = [obj.v_b_k_1; obj.s_dot_k_1];

            % Setting the parameter for the optimization solver
            obj.M = obj.casadi_optimizer.parameter(NDOF+6,NDOF+6);
            obj.h = obj.casadi_optimizer.parameter(NDOF+6);
            obj.Jc = obj.casadi_optimizer.parameter(12,NDOF+6);
            obj.Jcmm = obj.casadi_optimizer.parameter(6,NDOF+6);
            obj.nu_k = obj.casadi_optimizer.parameter(NDOF+6);
            obj.C_friction = obj.casadi_optimizer.parameter(38,12);
            obj.b_friction = obj.casadi_optimizer.parameter(38);
            obj.s_dot_star = obj.casadi_optimizer.parameter(NDOF);
            obj.H_star = obj.casadi_optimizer.parameter(6);
            obj.Jc_dot_nu = obj.casadi_optimizer.parameter(12);
            obj.tau_meas = obj.casadi_optimizer.parameter(NDOF);
            obj.f_meas = obj.casadi_optimizer.parameter(12);
            obj.s = obj.casadi_optimizer.parameter(NDOF);
            obj.CoM_measured = obj.casadi_optimizer.parameter(3);
            obj.K_friction = obj.casadi_optimizer.parameter(NDOF+6,NDOF+6);
            obj.t_step = obj.casadi_optimizer.parameter(1);
            obj.CoM_max = obj.casadi_optimizer.parameter(2);
            obj.CoM_min = obj.casadi_optimizer.parameter(2);
            obj.CoM_vel_max = obj.casadi_optimizer.parameter(2);
            obj.CoM_vel_min = obj.casadi_optimizer.parameter(2);
            obj.epsilon_CoM = obj.casadi_optimizer.parameter(2);
            obj.ActivateComHeightConstraint = obj.casadi_optimizer.parameter(1);
            obj.CoM_Z_reference_velocity = obj.casadi_optimizer.parameter(1);
            obj.feet_correction = obj.casadi_optimizer.parameter(12);
            obj.holonomicConstraintActivation = obj.casadi_optimizer.parameter(1);
            obj.JTorso= obj.casadi_optimizer.parameter(6,NDOF+6);
            obj.SaturationMaxJoints = obj.casadi_optimizer.parameter(NDOF);
            obj.SaturationMinJoints = obj.casadi_optimizer.parameter(NDOF);
            obj.J_com = obj.casadi_optimizer.parameter(6,NDOF+6);

            %Weigth
            Weigth.PosturalTask    = 100;
            Weigth.MomentumLinear  = 2.0;
            Weigth.MomentumAngular = 30.00;
            Weigth.RegTorques      = 0.001;
            Weigth.RegVelocities   = 128;
            Weigth.Wrenches        = 0.01;
            Weigth.HolonomicConstraint = 10;
            % Selector Matrix

            B           = [ zeros(6,NDOF);  ...
                            eye(NDOF,NDOF)];

            % Setting the objective function
            obj.casadi_optimizer.minimize(Weigth.RegTorques*sumsqr(obj.tau_k - obj.tau_meas)+             ... % Torque Regularization
                Weigth.RegVelocities*sumsqr(nu_k_1(1:6))+                                                 ... % Base Velocity Regularization
                Weigth.PosturalTask*sumsqr(obj.s_dot_k_1-obj.s_dot_star));                                ... % Postural Task
                %Weigth.HolonomicConstraint*(1-obj.holonomicConstraintActivation)*sumsqr(obj.Jc*(nu_k_1)));                                     ... %HolonomicCOnstraint
                %Weigth.Wrenches*sumsqr(obj.f_k) +                                                        ... %Wrenches
                %Weigth.MomentumAngular*sumsqr(obj.H_star(4:6)-obj.Jcmm(4:6,:)*nu_k_1));                  ... % Angular Momentum
                %Weigth.MomentumLinear*sumsqr(obj.H_star(1:3)-obj.Jcmm(1:3,:)*nu_k_1));                   ... % Linear Momentum

            % Setting Constraint
            % Dynamics
            obj.casadi_optimizer.subject_to((obj.M/obj.t_step + obj.K_friction)*(nu_k_1)-(obj.M*obj.nu_k)/obj.t_step+obj.h== B*obj.tau_k +obj.Jc'*obj.f_k);
            % Holonomic Constraint
            % acceleration-wise
            %obj.casadi_optimizer.subject_to(obj.holonomicConstraintActivation*(obj.Jc*(nu_k_1 - obj.nu_k)/obj.t_step) == obj.holonomicConstraintActivation*(-obj.Jc_dot_nu));
            % velocity-wise
            obj.casadi_optimizer.subject_to(obj.holonomicConstraintActivation*obj.Jc*(nu_k_1)==obj.holonomicConstraintActivation*obj.feet_correction);
            % Friction Cones
            obj.casadi_optimizer.subject_to(obj.C_friction*obj.f_k<obj.b_friction);

            %% Constraint on the CoM Velocity
            obj.casadi_optimizer.subject_to(obj.J_com(1,:)*nu_k_1<=(tanh(obj.epsilon_CoM*(obj.CoM_max(1)-obj.CoM_measured(1)))*obj.CoM_vel_max(1)));
            obj.casadi_optimizer.subject_to(obj.J_com(1,:)*nu_k_1>=(tanh(obj.epsilon_CoM*(obj.CoM_measured(1)-obj.CoM_min(1)))*obj.CoM_vel_min(1)));
            obj.casadi_optimizer.subject_to(obj.J_com(2,:)*nu_k_1<=(tanh(obj.epsilon_CoM*(obj.CoM_max(2)-obj.CoM_measured(2)))*obj.CoM_vel_max(2)));
            obj.casadi_optimizer.subject_to(obj.J_com(2,:)*nu_k_1>=(tanh(obj.epsilon_CoM*(obj.CoM_measured(2)-obj.CoM_min(2)))*obj.CoM_vel_min(2)));

            % Constraint the Torso velocity on the z axis
            %obj.casadi_optimizer.subject_to(obj.ActivateComHeightConstraint*obj.Jcmm(3,7:end)*nu_k_1(7:end) >= obj.ActivateComHeightConstraint*obj.M(1,1)*obj.CoM_Z_reference_velocity);
            obj.casadi_optimizer.subject_to(obj.ActivateComHeightConstraint*obj.JTorso(3,:)*nu_k_1 >= obj.ActivateComHeightConstraint*obj.CoM_Z_reference_velocity);
%
            % Angular Momentum
            % With Control Law
            % obj.casadi_optimizer.subject_to(obj.H_star(4:6,:)==obj.Jcmm(4:6,:)*nu_k_1);
            % Setting to Zero
            %obj.casadi_optimizer.subject_to(obj.Jcmm(4:6,:)*nu_k_1 == zeros(3,1));
            %THE BOUND FOR NOW ARE DE-ACTIVATED
            % Setting upper and lower bound
            % Lower and Upper Bound NOTE by using a casadi optimizer object, we loose the capability of interpreting the lower and
            % uppe bound in a smart way.
            % obj.casadi_optimizer.subject_to(Sat.BaseLinearVelocity_min<= obj.v_b_k_1(1:3)<=Sat.BaseLinearVelocity_max);
            % obj.casadi_optimizer.subject_to(Sat.BaseAngVelocity_min<= obj.v_b_k_1(4:end)<=Sat.BaseAngVelocity_max);
            %obj.casadi_optimizer.subject_to(obj.SaturationMinJoints<= obj.s_dot_k_1);
            %obj.casadi_optimizer.subject_to(obj.s_dot_k_1<= obj.SaturationMaxJoints);
            % obj.casadi_optimizeru.sbject_to(Sat.torques_min(1:3)<= obj.tau_k(1:3)<= Sat.torques_max(1:3));

        end

        function [v_b_star,s_dot_k_1_star,tau_star,f_star] = stepImpl(obj,M,h,Jc,Jcmm,nu_k,Jc_dot_nu,C_friction,b_friction,H_star,s_dot_star,tau_meas, contact_status, f_meas, x_com, J_com, w_H_b,...
                s,CoM_measured,ADD_FRICTION, K_viscous_friction, T, Gamma, baseVelocitySat, jointsVelocitySat, torqueSat, wrenchesSat, tStep, NDOF,epsilon_CoM, CoM_limits,CoM_vel_limits, activateCoMHeightTask, CoMZReferenceVelocity, correctionFeet, activateHolonomicConstraint, JTorso )
            solver_succeded = true;

            % Compute Centroidayl Dynamics
            selector = zeros(length(nu_k),1);
            selector(3) = 1 ;
            DYNAMICS.M = M;
            DYNAMICS.h = h;
            DYNAMICS.g = M*selector*(9.81);
            DYNAMICS.Jc = Jc;
            DYNAMICS.dJc_nu = Jc_dot_nu;
            FORKINEMATICS.xCoM = x_com;
            vel_com = J_com*nu_k;
            FORKINEMATICS.dxCoM = vel_com(1:3);

            STATE.nu = nu_k;
            STATE.dx_b = nu_k(1:3);
            STATE.x_b = w_H_b(1:3,4);
            centroidalDyn = centroidalConversion(DYNAMICS,FORKINEMATICS,STATE);

            %K friction
            K_friction_new = zeros(NDOF+6);

            if(ADD_FRICTION)
                invTGamma                = eye(size(Gamma))/(T*Gamma);
                invTGamma_t              = eye(size(Gamma))/(transpose(T*Gamma));
                K_friction_new(7:end, 7:end) = invTGamma_t*K_viscous_friction*invTGamma;
            end

            % Update casadi optimizer parameters
            obj.casadi_optimizer.set_value(obj.M, centroidalDyn.M);
            obj.casadi_optimizer.set_value(obj.h,( centroidalDyn.C_nu+centroidalDyn.g));
            obj.casadi_optimizer.set_value(obj.Jc, centroidalDyn.Jc);
            obj.casadi_optimizer.set_value(obj.Jcmm, Jcmm);
            obj.casadi_optimizer.set_value(obj.nu_k, centroidalDyn.nu);
            obj.casadi_optimizer.set_value(obj.C_friction,C_friction);
            obj.casadi_optimizer.set_value(obj.b_friction,b_friction);
            obj.casadi_optimizer.set_value(obj.s_dot_star,s_dot_star);
            obj.casadi_optimizer.set_value(obj.H_star,H_star);
            obj.casadi_optimizer.set_value(obj.Jc_dot_nu, centroidalDyn.dJc_nu);
            obj.casadi_optimizer.set_value(obj.tau_meas, tau_meas);
            obj.casadi_optimizer.set_value(obj.f_meas, f_meas);
            obj.casadi_optimizer.set_value(obj.s,s);
            obj.casadi_optimizer.set_value(obj.CoM_measured,CoM_measured)
            obj.casadi_optimizer.set_value(obj.K_friction,K_friction_new);
            obj.casadi_optimizer.set_value(obj.t_step,tStep);
            obj.casadi_optimizer.set_value(obj.CoM_max, CoM_limits(:,2));
            obj.casadi_optimizer.set_value(obj.CoM_min, CoM_limits(:,1));
            obj.casadi_optimizer.set_value(obj.CoM_vel_max, CoM_vel_limits(:,2));
            obj.casadi_optimizer.set_value(obj.CoM_vel_min, CoM_vel_limits(:,1));
            obj.casadi_optimizer.set_value(obj.epsilon_CoM, epsilon_CoM);
            obj.casadi_optimizer.set_value(obj.ActivateComHeightConstraint,activateCoMHeightTask);
            obj.casadi_optimizer.set_value(obj.CoM_Z_reference_velocity, CoMZReferenceVelocity);
            obj.casadi_optimizer.set_value(obj.feet_correction, correctionFeet);
            obj.casadi_optimizer.set_value(obj.holonomicConstraintActivation, activateHolonomicConstraint);
            obj.casadi_optimizer.set_value(obj.JTorso,JTorso);
            obj.casadi_optimizer.set_value(obj.J_com,J_com);
%             obj.casadi_optimizer.set_value(obj.SaturationMaxJoints,jointsVelocitySat(:,2));
%             obj.casadi_optimizer.set_value(obj.SaturationMinJoints, jointsVelocitySat(:,1));

            % Computing the solution
            try
                obj.sol = obj.casadi_optimizer.solve; % could raise an error
                v_b_star = full(obj.sol.value(obj.v_b_k_1));
                s_dot_k_1_star = full(obj.sol.value(obj.s_dot_k_1));
                tau_star = full(obj.sol.value(obj.tau_k));
                f_star = full(obj.sol.value(obj.f_k));
                obj.casadi_optimizer.set_initial(obj.casadi_optimizer.x, obj.sol.value(obj.casadi_optimizer.x));
            catch exception
                try
                    solver_succeded = false;
                    obj.casadi_optimizer.debug.show_infeasibilities;
                    disp(exception.getReport)
                catch
                    disp(exception.message)
                end
            end

            % If the solver did not succed, put to zero the control input
            % and increase the related counter
            if(~solver_succeded)
                obj.solver_fails_counter = obj.solver_fails_counter +1;
                v_b_star = zeros(6,1);
                s_dot_k_1_star = zeros(26,1);
                tau_star = zeros(26,1);
                f_star = zeros(12,1);
                if(obj.solver_fails_counter>6)
                    error('solver failed more that 4 times');
                end
            else
                obj.solver_fails_counter = 0;

            end

        end

        function resetImpl(obj)
            % Initialize discrete-state properties.
        end
    end
end

author: robotology

repo: whole-body-controllers


text:
function [noFluxRxns, noFluxRxnsRelaxed, subGraphs, notProducedMets, minToConnect,...
    neededForProductionMat, canProduceWithoutInput, canConsumeWithoutOutput, ...
    connectedFromTemplates, addedFromTemplates]=gapReport(model, templateModels)
% gapReport
%   Performs a gap analysis and summarizes the results
%
%   model                       a model structure
%   templateModels              a cell array of template models to use for
%                               gap filling (opt)
%
%   noFluxRxns                  cell array with reactions that cannot carry
%                               flux
%   noFluxRxnsRelaxed           cell array with reactions that cannot carry
%                               flux even if the mass balance constraint is
%                               relaxed so that it is allowed to have
%                               net production of all metabolites
%   subGraphs                   structure with the metabolites in each of
%                               the isolated sub networks
%   notProducedMets             cell array with the metabolites that
%                               couldn't have net production
%   minToConnect                structure with the minimal number of
%                               metabolites that need to be connected in
%                               order to be able to produce all other
%                               metabolites and which metabolites each of
%                               them connects
%   neededForProductionMat      matrix where n x m is true if metabolite n
%                               allows for production of metabolite m
%   canProduceWithoutInput      cell array with metabolites that could be
%                               produced even when there is no input to the
%                               model
%   canConsumeWithoutOutput     cell array with metabolites that could be
%                               consumed even when there is no output from
%                               the model
%   connectedFromTemplates      cell array with the reactions that could be
%                               connected using the template models
%   addedFromTemplates          structure with the reactions that were
%                               added from the template models and which
%                               model they were added from
%
%   Usage: [noFluxRxns, noFluxRxnsRelaxed, subGraphs, notProducedMets, minToConnect,...
%    neededForProductionMat, connectedFromTemplates, addedFromTemplates]=...
%    gapReport(model, templateModels)

if nargin<2
    templateModels=[];
    connectedFromTemplates=[];
    addedFromTemplates=[];
end

fprintf(['Gap analysis for ' model.id ' - ' model.name '\n\n']);
if isfield(model,'unconstrained')
    calculateINOUT=true;
    closedModel=model;
    model=simplifyModel(model);
else
    canConsumeWithoutOutput={};
    canProduceWithoutInput={};
    calculateINOUT=false;
end

model2=model;
model2.b=[model2.b inf(numel(model2.mets),1)];
I=haveFlux(model);
noFluxRxns=model.rxns(~I);
J=haveFlux(model2);
noFluxRxnsRelaxed=model2.rxns(~J);
bModel=removeReactions(model,~I,true,true);
cModel=removeReactions(model2,~J,true,true);
fprintf('***Overview\n');
fprintf([num2str(numel(model.rxns)-sum(I)) ' out of ' num2str(numel(model.rxns))...
    ' reactions cannot carry flux (' num2str(numel(model.rxns)-sum(J)) ' if net production of all metabolites is allowed)\n']);
fprintf([num2str(numel(model.mets)-numel(bModel.mets)) ' out of ' num2str(numel(model.mets))...
    ' metabolites are unreachable (' num2str(numel(model.mets)-numel(cModel.mets)) ' if net production of all metabolites is allowed)\n']);

fprintf('\n***Isolated subnetworks\n');
subGraphs=getAllSubGraphs(model);
fprintf(['A total of ' num2str(size(subGraphs,2)) ' isolated sub-networks are present in the model\n']);
for i=1:size(subGraphs,2)
    fprintf(['\t' num2str(i) '. ' num2str(sum(subGraphs(:,i))) ' metabolites\n']);
end

fprintf('\n***Metabolite connectivity\n');
[notProducedMets, ~, neededForProductionMat,minToConnect]=checkProduction(model,true,model.comps,false);
fprintf(['To enable net production of all metabolites, a total of ' num2str(numel(minToConnect)) ' metabolites must be connected\n']);
fprintf('Top 10 metabolites to connect:\n');
for i=1:min(10,numel(minToConnect))
    fprintf(['\t' num2str(i) '. ' minToConnect{i} '\n']);
end

if calculateINOUT==true
    fprintf('\n***Mass balancing\n');
    produced=canProduce(closedModel);
    canProduceWithoutInput=closedModel.mets(produced);
    consumed=canConsume(closedModel);
    canConsumeWithoutOutput=closedModel.mets(consumed);
    fprintf([num2str(numel(canConsumeWithoutOutput)) ' metabolites could be consumed without any outputs\n' num2str(numel(canProduceWithoutInput)) ' metabolites could be produced without any inputs\n']);
end

if ~isempty(templateModels)
    fprintf('\n***Automated gap-filling\n');
    [connectedFromTemplates, ~, addedFromTemplates]=fillGaps(model,templateModels);
    t=templateModels{1}.id;
    for i=2:numel(templateModels)
        t=[t ', ' templateModels{i}.id];
    end
    fprintf([num2str(numel(connectedFromTemplates)) ' unconnected reactions can be connected by including ' num2str(numel(addedFromTemplates)) ' reactions from\n' t '\n']);
end
end

author: SysBioChalmers

repo: RAVEN