text:
function Fog_simulation_Cityscapes(task_id, dataset_split, refinement_level,...
    variant, beta, output_root_directory, images_per_task)
%FOG_SIMULATION_CITYSCAPES  Simulate fog for a batch of images from Cityscapes
%and write the results. Structured for execution on a cluster.
%
%   INPUTS:
%
%   -|task_id|: ID of the task. Used in conjuction with |images_per_task| to
%    determine which images out of the entire dataset will form the batch that
%    will be processed by this task.
%
%   -|dataset_split|: string that indicates which subset of Cityscapes is used
%    for running the simulation, e.g. 'trainval', 'train_extra'.
%
%   -|refinement_level|: string that indicates the refinement that is applied to
%    the used subset, e.g. 'full', 'refined'.
%
%   -|variant|: string that indicates the type of fog simulation that is being
%    run, e.g. 'stereoscopic_inpainting_with_guided_filtering',
%    'nearest_neighbor'.
%
%   -|beta|: value of attenuation coefficient that determines fog density.
%
%   -|cityscapes_output_root_directory|: full path to directory under which the
%    results of the simulation are written, e.g.
%    '/scratch_net/nowin/csakarid/Toyota-foggy/data/Cityscapes'.
%
%   -|images_per_task|: maximum number of images for which fog simulation is run
%    in each task.
%    19997 train_extra Cityscapes images -> 200 images per task * 100 tasks.
%    5000 trainvaltest Cityscapes images -> 100 images per task * 50 tasks.
%    550 trainval refined Cityscapes images -> 22 images per task * 25 tasks.

if ischar(task_id)
    task_id = str2double(task_id);
end

% ------------------------------------------------------------------------------

% Check which variant of the simulation is being run and distinguish between
% different cases.

% Get abbreviations to be included in output names.
abbreviations_of_variants =...
    containers.Map({'stereoscopic_inpainting_with_guided_filtering',...
    'nearest_neighbor'}, {'stereogf', 'nearest'});
variant_abbreviation = abbreviations_of_variants(variant);

% ------------------------------------------------------------------------------

% Add paths to functions that are called for fog simulation.

current_script_full_name = mfilename('fullpath');
current_script_directory = fileparts(current_script_full_name);
addpath(fullfile(current_script_directory, '..', 'utilities'));
addpath_relative_to_caller(current_script_full_name,...
    fullfile('..', 'Fog_simulation'));
addpath_relative_to_caller(current_script_full_name,...
    fullfile('..', 'Depth_processing'));
addpath_relative_to_caller(current_script_full_name,...
    fullfile('..', 'Dehazing', 'Dark_channel_prior'));
addpath_relative_to_caller(current_script_full_name,...
    fullfile('..', 'external', 'SLIC_mex'));

% ------------------------------------------------------------------------------

% Create lists of input files for various modalities.

cityscapes_file_lists_directory = fullfile(current_script_directory, '..',...
    '..', 'data', 'lists_file_names');

% Left input images.
list_of_left_images_file = strcat('leftImg8bit_', dataset_split, '_',...
    refinement_level, '_clean_filenames.txt');
fid = fopen(fullfile(cityscapes_file_lists_directory,...
    list_of_left_images_file));
left_image_file_names = textscan(fid, '%s');
fclose(fid);
left_image_file_names = left_image_file_names{1};

number_of_images = length(left_image_file_names);

% Right input images.
list_of_right_images_file = strcat('rightImg8bit_', dataset_split, '_',...
    refinement_level, '_filenames.txt');
fid = fopen(fullfile(cityscapes_file_lists_directory,...
    list_of_right_images_file));
right_image_file_names = textscan(fid, '%s');
fclose(fid);
right_image_file_names = right_image_file_names{1};

% Disparity map files.
list_of_disparity_maps_file = strcat('disparity_', dataset_split, '_',...
    refinement_level, '_filenames.txt');
fid = fopen(fullfile(cityscapes_file_lists_directory,...
    list_of_disparity_maps_file));
disparity_file_names = textscan(fid, '%s');
fclose(fid);
disparity_file_names = disparity_file_names{1};

% Camera parameters files.
list_of_camera_parameters_file = strcat('camera_', dataset_split, '_',...
    refinement_level, '_filenames.txt');
fid = fopen(fullfile(cityscapes_file_lists_directory,...
    list_of_camera_parameters_file));
camera_parameters_file_names = textscan(fid, '%s');
fclose(fid);
camera_parameters_file_names = camera_parameters_file_names{1};

% Based on the selected simulation variant, form the name of the corresponding
% text file with the file names of the respective depth maps.
switch variant_abbreviation
    case 'stereogf'
        list_of_depth_maps_file = strcat('depth_stereoscopic_',...
            dataset_split, '_', refinement_level, '_filenames.txt');
    case 'nearest'
        list_of_depth_maps_file = strcat('depth_nearest_', dataset_split,...
            '_', refinement_level, '_filenames.txt');
end

% Get all names of depth map files, whether they already exist or not.
fid = fopen(fullfile(cityscapes_file_lists_directory,...
    list_of_depth_maps_file));
depth_file_names = textscan(fid, '%s');
fclose(fid);
depth_file_names = depth_file_names{1};

% Sanity checks. Total number of files should be identical for all modalities.
assert(number_of_images == length(right_image_file_names));
assert(number_of_images == length(disparity_file_names));
assert(number_of_images == length(camera_parameters_file_names));
assert(number_of_images == length(depth_file_names));

% ------------------------------------------------------------------------------

% Determine current batch.

% Determine the set of images that are to be processed in the current task.
batch_ind = (task_id - 1) * images_per_task + 1:task_id * images_per_task;
if batch_ind(1) > number_of_images
    return;
end
if batch_ind(end) > number_of_images
    % Truncate for last task.
    batch_ind = batch_ind(1:number_of_images - batch_ind(1) + 1);
end

% ------------------------------------------------------------------------------

% Instantiate individual components of the fog simulation pipeline that are
% required for the selected variant and set their parameters.

% Attenuation coefficient. Set it to meaningful values.
attenuation_coefficient_method = @scattering_coefficient_fixed;
beta_parameters.beta = beta;

% Atmospheric light.
atmospheric_light_estimation = @estimate_atmospheric_light_rf;

% Select the optical model which is used for the synthesis of foggy images from
% the clean input images.
fog_optical_model = @haze_linear;

% Transmittance model of the clear scene radiance.
transmittance_model = @transmission_homogeneous_medium;

% Specify whether depth maps will be loaded from existing .mat files if possible
% or computed on the fly.
load_depth = true;

% Specify whether the depth maps used for fog simulation will be saved to .mat
% files.
save_depth = true;

% Select the method for depth map completion and denoising, and the method for
% subsequent postprocessing of the initial transmittance map.
switch variant_abbreviation
    case 'stereogf'
        depth_completion_method =...
            @depth_in_meters_cityscapes_stereoscopic_inpainting;
        transmittance_postprocessing =...
            @transmission_postprocessing_guided_filter;
    case 'nearest'
        depth_completion_method = @depth_in_meters_cityscapes_nn_inpainting;
end

% ------------------------------------------------------------------------------

% Output specifications.

% Specify output format for synthetic foggy images and corresponding
% transmittance maps.
output_format = '.png';

% Determine initial part of the path of output files.
cityscapes_output_directories_basename = strcat('leftImg8bit_',...
    dataset_split, '_', refinement_level, '_', variant_abbreviation,...
    '_beta_', num2str(beta_parameters.beta));
cityscapes_output_foggy_directory = fullfile(output_root_directory,...
    strcat(cityscapes_output_directories_basename, '_foggy'));
cityscapes_output_transmittance_directory = fullfile(output_root_directory,...
    strcat(cityscapes_output_directories_basename, '_transmittance'));

% ------------------------------------------------------------------------------

% Run fog simulation on current batch using the specified settings.

switch variant_abbreviation

    % Stereoscopic inpainting followed by guided filtering.
    case 'stereogf'
        clean2foggy_Cityscapes(left_image_file_names(batch_ind),...
            right_image_file_names(batch_ind),...
            disparity_file_names(batch_ind),...
            camera_parameters_file_names(batch_ind),...
            depth_file_names(batch_ind), load_depth, save_depth,...
            depth_completion_method, attenuation_coefficient_method,...
            beta_parameters, atmospheric_light_estimation,...
            transmittance_model, transmittance_postprocessing,...
            fog_optical_model, output_root_directory,...
            cityscapes_output_foggy_directory,...
            cityscapes_output_transmittance_directory, output_format);

    % Nearest neighbor inpainting of depth values.
    case 'nearest'
        clean2foggy_Cityscapes_nn(left_image_file_names(batch_ind),...
            disparity_file_names(batch_ind),...
            camera_parameters_file_names(batch_ind),...
            attenuation_coefficient_method, beta_parameters,...
            atmospheric_light_estimation, transmittance_model,...
            fog_optical_model, cityscapes_output_foggy_directory,...
            cityscapes_output_transmittance_directory, output_format);
end

end

author: sakaridis


text:
function [freq] = ft_appendfreq(cfg, varargin)

% FT_APPENDFREQ concatenates multiple frequency or time-frequency data structures
% that have been processed separately. If the input data structures contain different
% channels, it will be concatenated along the channel direction. If the channels are
% identical in the input data structures, the data will be concatenated along the
% repetition dimension.
%
% Use as
%  combined = ft_appendfreq(cfg, freq1, freq2, ...)
%
% The configuration should contain
%   cfg.parameter  = string, the name of the field to concatenate
%
% The configuration can optionally contain
%   cfg.appenddim  = string, the dimension to concatenate over (default is automatic)
%   cfg.tolerance  = scalar, tolerance to determine how different the frequency and/or
%                    time axes are allowed to still be considered compatible (default = 1e-5)
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a
% *.mat file on disk and/or the output data will be written to a *.mat file.
% These mat files should contain only a single variable, corresponding with
% the input/output structure.
%
% See also FT_FREQANALYSIS, FT_DATATYPE_FREQ, FT_APPENDDATA, FT_APPENDTIMELOCK,
% FT_APPENDSENS

% Copyright (C) 2011-2017, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar    varargin
ft_preamble provenance varargin

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
for i=1:length(varargin)
  varargin{i} = ft_checkdata(varargin{i}, 'datatype', {'freq', 'freq+comp'}, 'feedback', 'yes');
end

% set the defaults
cfg.channel    = ft_getopt(cfg, 'channel', 'all');
cfg.channelcmb = ft_getopt(cfg, 'channelcmb', {'all' 'all'});
cfg.parameter  = ft_getopt(cfg, 'parameter', []);
cfg.appenddim  = ft_getopt(cfg, 'appenddim', []);
cfg.tolerance  = ft_getopt(cfg, 'tolerance',  1e-5); % this is passed to append_common, which passes it to ft_selectdata
cfg.appendsens = ft_getopt(cfg, 'appendsens', 'no');

hastime = isfield(varargin{1}, 'time');
hasfreq = isfield(varargin{1}, 'freq');

if isempty(cfg.appenddim) || strcmp(cfg.appenddim, 'auto')
  if hastime && hasfreq
    if checkchan(varargin{:}, 'identical') && checkfreq(varargin{:}, 'identical', cfg.tolerance) && checktime(varargin{:}, 'identical', cfg.tolerance)
      cfg.appenddim = 'rpt';
    elseif checkfreq(varargin{:}, 'unique', cfg.tolerance) && checktime(varargin{:}, 'identical', cfg.tolerance)
      cfg.appenddim = 'freq';
    elseif checktime(varargin{:}, 'unique', cfg.tolerance) && checkfreq(varargin{:}, 'identical', cfg.tolerance)
      cfg.appenddim = 'time';
    elseif checkchan(varargin{:}, 'unique')
      cfg.appenddim = 'chan';
    else
      ft_error('cfg.appenddim should be specified');
    end
  else
    if checkchan(varargin{:}, 'identical') && checkfreq(varargin{:}, 'identical', cfg.tolerance)
      cfg.appenddim = 'rpt';
    elseif checkfreq(varargin{:}, 'unique', cfg.tolerance)
      cfg.appenddim = 'freq';
    elseif checkchan(varargin{:}, 'unique')
      cfg.appenddim = 'chan';
    else
      ft_error('cfg.appenddim should be specified');
    end
  end
end
ft_info('concatenating over the "%s" dimension\n', cfg.appenddim);

if isempty(cfg.parameter)
  fn = fieldnames(varargin{1});
  for i=2:numel(varargin)
    fn = intersect(fn, fieldnames(varargin{i}));
  end
  cfg.parameter = setdiff(fn, ignorefields('appendfreq'));
elseif ischar(cfg.parameter)
  cfg.parameter = {cfg.parameter};
end
assert(~isempty(cfg.parameter), 'cfg.parameter should be specified');

% use a low-level function that is shared with the other ft_appendxxx functions
freq = append_common(cfg, varargin{:});

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   varargin
ft_postamble provenance freq
ft_postamble history    freq
ft_postamble savevar    freq

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function polygon_properties_test05 ( )

%*****************************************************************************80
%
%% POLYGON_PROPERTIES_TEST05 tests POLYGON_CONTAINS_POINT, POLYGON_CONTAINS_POINT_2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2014
%
%  Author:
%
%    John Burkardt
%
  n = 5;
  test_num = 4;
  p_test = [ ...
    1.0,  1.0; ...
    3.0,  4.0; ...
    0.0,  2.0; ...
    0.5, -0.25 ]';
  v = [ ...
    0.0, 0.0; ...
    1.0, 0.0; ...
    2.0, 1.0; ...
    1.0, 2.0; ...
    0.0, 2.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'POLYGON_PROPERTIES_TEST05\n' );
  fprintf ( 1, '  POLYGON_CONTAINS_POINT determines if\n' );
  fprintf ( 1, '  a point is in a polygon.\n' );
  fprintf ( 1, '  POLYGON_CONTAINS_POINT_2 determines if\n' );
  fprintf ( 1, '  a point is in a polygon.\n' );

  r8mat_transpose_print ( 2, n, v, '  The polygon vertices:' );

  fprintf ( 1, '\n' );
  fprintf ( 1, '          P          In1  In2\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num

    p(1:2,1) = p_test(1:2,test);

    inside1 = polygon_contains_point ( n, v, p );

    inside2 = polygon_contains_point_2 ( n, v, p );

    fprintf ( 1, '  %14.6g%14.6g    %d    %d\n', p(1:2), inside1, inside2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = laplace_pdf(x,mu,sigma)
%LAPLACE_PDF     Laplace probability density function (pdf).
%
%   Y = LAPLACE_PDF(X,MU,SIGMA) Returns the Laplace pdf with
%   mean, MU, and scale, SIGMA, at the values in X.
%
%   The size of Y is the common size of the input arguments. A scalar input
%   functions as a constant matrix of the same size as the other inputs.
%
%   Default values for MU and SIGMA are 0 and 1 respectively.

% Copyright (c) 2005 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 3,
  sigma = 1;
end

if nargin < 2;
  mu = 0;
end

if nargin < 1,
  error('Requires at least one input argument.');
end

y= -abs(x-mu)./sigma - log(2.*sigma);
y=exp(y);


text:
function rates
    % NO DESCRIPTION
    % 1: without declustering, 2: with declustering
    report_this_filefun()

    cd ~/zmap/src/thomas/seismicrates

    N=1
    for i=1:N
        params=sr_startZ(1);
    end
end

author: CelsoReyes


text:
function img=ndimfilter(im,kernel,varargin)
%
% img=ndimfilter(im,kernel,r,sigma)
%
% filter an ND array using a specified filter using convolution
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%    im: input ND array
%    kernel: can be an ND array, or a string. if string, the below filters
%        are supported:
%        'box': box filter (need r)
%        'gaussian': Gaussian filter (need r,sigma input)
%    r: kernel half-width, the output is 2*r+1 in each dimension; if
%       missing, use 1
%    sigma: the standard deviation of the Gaussian; if not given, use 1; if
%       set to inf, output box filter
%
% output:
%    img: the filtered ND array
%
% -- this function is part of the Iso2Mesh Toolbox (http://iso2mesh.sf.net)
%    License: GPL v3 or later, see LICENSE.txt for details
%

if(nargin<2)
    kernel='box';
end
if(ischar(kernel))
    switch(kernel)
        case 'box'
            kernel=ndgaussian(varargin{1},inf,ndims(im));
        case 'gaussian'
            kernel=ndgaussian(varargin{1},varargin{2},ndims(im));
        otherwise
            error('filter type %s is not supported',type);
    end
end

img=convn(im,kernel,'same');

author: fangq


text:
function inside = triangle_contains_point_2d_1 ( t, p )

%*****************************************************************************80
%
%% TRIANGLE_CONTAINS_POINT_2D_1 finds if a point is inside a triangle in 2D.
%
%  Discussion:
%
%    It is conventional to list the triangle vertices in counter clockwise
%    order.  However, this routine does not require a particular order
%    for the vertices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, logical INSIDE, is TRUE if the point is inside
%    the triangle or on its boundary.
%
  xsi = triangle_barycentric_2d ( t, p );

  if ( any ( xsi(1:3,1) < 0.0 ) )
    inside = 0;
  else
    inside = 1;
  end

  return
end


text:
function [manifest] = write_manifest(tracker, directory)
% write_manifest Write a manifest file for the tracker
%
% Write a manifest file for the tracker. The manifest file contains some system
% statistics that can be used for debugging and analysis.
%
% Input:
% - tracker: Tracker structure.
% - tracker: Optional directory (default is tracker result root).
%
% Output:
% - manifest: Full path to the manifest file that is located in tracker result directory.

if nargin < 2
    directory = tracker.directory;
end;

mkpath(directory);
manifest = fullfile(directory, 'manifest.txt');

[platform_str, platform_maxsize, platform_endian] = computer();

if is_octave()
    environment = 'octave';
else
    environment = 'matlab';
end;

environment_version = version();

fid = fopen(manifest, 'w');

votversion = toolkit_version();

fprintf(fid, 'toolkit.version=%d.%d.%d\n', votversion.major, votversion.minor, votversion.patch);
if isfield(votversion, 'build')
fprintf(fid, 'toolkit.build=%d\n', votversion.build);
end;
fprintf(fid, 'tracker.identifier=%s\n', tracker.identifier);
fprintf(fid, 'tracker.interpreter=%s\n', tracker.interpreter);
fprintf(fid, 'timestamp=%s\n', datestr(now, 31));
fprintf(fid, 'platform=%s\n', platform_str);
fprintf(fid, 'platform.maxsize=%f\n', platform_maxsize);
fprintf(fid, 'platform.endian=%s\n', platform_endian);
fprintf(fid, 'environment=%s\n', environment);
fprintf(fid, 'environment.version=%s\n', environment_version);

fclose(fid);

author: votchallenge

repo: toolkit-legacy


text:
function [par] = readprocpar(fpath, ext)

% READPROCPAR reads in parameter used in FID acquisition
%
%   PAR = READPROCPAR(FPATH) reads the parameters in the the
%   procpar file in the FID directory FPATH
%
% see also WRITEPROCPAR, READFID, WRITEFID

% Written by L Martyn Klassen
% Copyright 2003 Robarts Research Institute

if nargin < 1
   error('CFMM:readprocpar','One input argument required.');
end

if nargin < 2
   ext = [];
end

[data,fpath] = readprocpar2(fpath, ext);

% Initialize outputs
par.path     = fpath;
par.dp       = [];
par.seqfil   = [];
par.nf       = [];
par.ni       = [];
par.sw       = [];
par.ns       = [];
par.ne1      = [];
par.ne2      = [];
par.rfspoil  = [];
par.rfphase  = [];
par.nt       = [];
par.np       = [];
par.nv       = [];
par.nv2      = [];
par.nv3      = [];
par.ne       = [];
par.polarity = [];
par.evenecho = [];
par.tr       = [];
par.te       = [];
par.esp      = [];
par.espincr  = [];
par.nprof    = [];
par.tproj    = [];
par.phi      = [];
par.psi      = [];
par.theta    = [];
par.vpsi     = [];
par.vphi     = [];
par.vtheta   = [];
par.array{1} = [];
par.arraydim = [];
par.seqcon   = [];
par.lro      = [];
par.lpe      = [];
par.lpe2     = [];
par.pro      = [];
par.ppe      = [];
par.ppe2     = [];
par.pss      = [];
par.thk      = [];
par.thk2     = [];
par.pos1     = [];
par.pos2     = [];
par.pos3     = [];
par.vox1     = [];
par.vox2     = [];
par.vox3     = [];
par.nD       = [];
par.cntr     = [];
par.gain     = [];
par.shimset  = [];
par.z1       = [];
par.z2       = [];
par.z3       = [];
par.z4       = [];
par.z5       = [];
par.z6       = [];
par.z7       = [];
par.z8       = [];
par.x1       = [];
par.y1       = [];
par.xz       = [];
par.yz       = [];
par.xy       = [];
par.x3       = [];
par.y3       = [];
par.x4       = [];
par.y4       = [];
par.z1c      = [];
par.z2c      = [];
par.z3c      = [];
par.z4c      = [];
par.xz2      = [];
par.yz2      = [];
par.xz2      = [];
par.yz2      = [];
par.zxy      = [];
par.z3x      = [];
par.z3y      = [];
par.zx3      = [];
par.zy3      = [];
par.z4x        = [];
par.z4y        = [];
par.z5x        = [];
par.z5y        = [];
par.x2y2       = [];
par.z2xy       = [];
par.z3xy       = [];
par.z2x3       = [];
par.z2y3       = [];
par.z3x3       = [];
par.z3y3       = [];
par.z4xy       = [];
par.zx2y2      = [];
par.z2x2y2     = [];
par.z3x2y2     = [];
par.z4x2y2     = [];
par.petable    = [];
par.nrcvrs     = [];
par.trise      = [];
par.at         = [];
par.gro        = [];
par.gmax       = [];
par.intlv      = [];
par.rcvrs      = [];
par.celem      = [];
par.arrayelemts = [];
par.contrast   = [];
par.tep        = [];
par.date       = [];
par.ti         = [];
par.gss2       = [];
par.gss        = [];
par.tpwri      = [];
par.tpwr1      = [];
par.tpwr2      = [];
par.orient     = [];
par.tof        = [];
par.resto      = [];
par.grox       = [];
par.groy       = [];
par.fov        = [];
par.res        = [];
par.npix       = [];
par.nseg       = [];
par.nzseg      = [];
par.waveform   = [];
par.SR         = [];
par.gradfrac   = [];
par.sfrq       = [];
par.B0         = [];
par.dtmap      = [];
par.nnav       = [];
par.tnav       = [];
par.fast       = [];
par.bt         = [];
par.nhomo      = [];
par.fpmult     = [];
par.d1         = [];
par.ss         = [];
par.ssc        = [];
par.r1         = [];
par.r2         = [];
par.ps_coils   = [];
par.coil_array = [];
par.nav        = [];
par.fliplist   = [];
par.varflip    = [];
par.nfreq      = [];
par.freq       = [];
par.flip       = [];
par.flip1      = [];
par.flipprep   = [];
par.seg        = [];
par.state      = [];
par.rfdelay    = [];
par.gro        = [];
par.gimp       = [];
par.SR         = [];
par.readaxis   = [];
par.timescale  = [];
par.etl        = [];
par.grof       = [];
par.Po         = [];
par.Psl        = [];
par.console    = [];
par.shimcoils  = [];
par.spiral_gmax      = [];
par.spiral_gamma     = [];
par.spiral_delay     = [];
par.spiral_tep       = [];
par.dtg              = [];
par.nturns           = [];
par.direction        = [];
par.ninterleave      = [];
par.tn               = [];
par.randomseed       = [];
par.interleave_order = [];
par.profile          = [];
par.image            = [];
par.spiral_version   = [];
par.spiral_filter    = [];
par.spiral_alpha     = [];
par.spiral_density   = [];
par.navigator        = [];
par.tfirst           = [];
par.tpe              = [];
par.ky_order         = [];
par.alternate        = [];
par.offlineAverages  = [];
par.threshold        = [];
par.cluster          = [];
par.weightfit        = [];
par.offlineAverage   = [];
par.bipolar          = [];
par.tpwr             = [];
par.tpwrf            = [];
par.dpwr             = [];
par.dpwrf            = [];
par.echoes           = [];
par.combinedAcq      = [];
par.reductionFactor  = [];
par.referenceUnits   = [];
par.referenceLines   = [];

names = fieldnames(par);
params = fieldnames(data);

for i = 1:length(params);
   par.(params{i}) = data.(params{i}).value;
end

% Make the returned values fully backwards compatible.
for i = 1:length(names)
   if isfield(data,names{i})
      field = data.(names{i});
      switch field.basictype
         case 1
            par.(names{i}) = field.value;
         case {0 2}
            par.(names{i}) = field.value{1};
      end
   end
end

if isfield(data, 'rcvrs')
   par.nrcvrs = data.rcvrs.number;
end

if isfield(data, 'array')
   par.array = data.array.loops;
end

if isfield(data, 'shimcoils')
   par.shimcoils = data.shimcoils.value;
end

if isfield(data, 'ps_coils')
   par.ps_coils = data.ps_coils.value;
end

if isfield(data, 'coil_array')
   par.coil_array = data.coil_array.value;
end

par.path = fpath;

return

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function [W,H] = nmf_cjlin(V,Winit,Hinit,tol,timelimit,maxiter)
%
% Copyright (c) 2005-2006 Chih-Jen Lin
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions
% are met:
%
% 1. Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in the
% documentation and/or other materials provided with the distribution.
%
% 3. Neither name of copyright holders nor the names of its contributors
% may be used to endorse or promote products derived from this software
% without specific prior written permission.
%
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
% ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
% A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
%
% NMF by alternative non-negative least squares using projected gradients
% Author: Chih-Jen Lin, National Taiwan University
%
% W,H: output solution
% Winit,Hinit: initial solution
% tol: tolerance for a relative stopping condition
% timelimit, maxiter: limit of time and iterations

W = Winit; H = Hinit; initt = cputime;

gradW = W*(H*H') - V*H'; gradH = (W'*W)*H - W'*V;
initgrad = norm([gradW; gradH'],'fro');
fprintf('Init gradient norm %f\n', initgrad);
tolW = max(0.001,tol)*initgrad; tolH = tolW;

for iter=1:maxiter,
  % stopping condition
  projnorm = norm([gradW(gradW<0 | W>0); gradH(gradH<0 | H>0)]);
  if projnorm < tol*initgrad | cputime-initt > timelimit,
    break;
  end

  [W,gradW,iterW] = nlssubprob(V',H',W',tolW,1000); W = W'; gradW = gradW';
  if iterW==1,
    tolW = 0.1 * tolW;
  end

  [H,gradH,iterH] = nlssubprob(V,W,H,tolH,1000);
  if iterH==1,
    tolH = 0.1 * tolH;
  end

  if (iterW==1 & iterH==1 & tolH + tolW < tol*initgrad),
    fprintf('Failed to move\n'); break;
  end
  if rem(iter,10)==0, fprintf('.'); end
end
fprintf('\nIter = %d Final proj-grad norm %f\n', iter, projnorm);

function [H,grad,iter] = nlssubprob(V,W,Hinit,tol,maxiter)
% H, grad: output solution and gradient
% iter: #iterations used
% V, W: constant matrices
% Hinit: initial solution
% tol: stopping tolerance
% maxiter: limit of iterations

H = Hinit;
WtV = W'*V;
WtW = W'*W;

alpha = 1; beta = 0.1;
for iter=1:maxiter,
  grad = WtW*H - WtV;
  projgrad = norm(grad(grad < 0 | H >0));
  if projgrad < tol,
    break
  end

  % search step size
  Hn = max(H - alpha*grad, 0); d = Hn-H;
  gradd=sum(sum(grad.*d)); dQd = sum(sum((WtW*d).*d));
  if gradd + 0.5*dQd > 0.01*gradd,
    % decrease alpha
    while 1,
      alpha = alpha*beta;
      Hn = max(H - alpha*grad, 0); d = Hn-H;
      gradd=sum(sum(grad.*d)); dQd = sum(sum((WtW*d).*d));
      if gradd + 0.5*dQd <= 0.01*gradd | alpha < 1e-20,
        H = Hn; break;
      end
    end
  else
    % increase alpha
    while 1,
      Hp = Hn;
      alpha = alpha/beta;
      Hn = max(H - alpha*grad, 0); d = Hn-H;
      gradd=sum(sum(grad.*d)); dQd = sum(sum((WtW*d).*d));
      if gradd + 0.5*dQd > 0.01*gradd | Hn == Hp | alpha > 1e10,
        H = Hp; alpha = alpha*beta; break;
      end
    end
  end
end

if iter==maxiter,
  fprintf('Max iter in nlssubprob\n');
end

author: andrewssobral

repo: lrslibrary


text:
% Gianni Schena  July 2005, schena@units.it
% Lattice Boltzmann LBE, geometry: D2Q9, model: BGK
% Application to permeability in porous media (low porosity)

Restart=false
logical(Restart);

if Restart==false;
close all, clear all % start from scratch and clean ...
Restart=false;
Pois_test=false, % for a test without obstacles
tic
%   IN
% |vvvv|    + y
% |vvvv|     ^
% |vvvv|     | -> + x
%  OUT

% Pores in 2D : Wet and Dry locations (Wet ==1 , Dry ==0 )
wXh_Dry=[3,1];wXh_Wet=[3,5];
%wXh_Dry=[3,0];wXh_Wet=[3,5]; % test non obstacles !!!
A=repmat([zeros(wXh_Dry),ones(wXh_Wet)],[1,3]);A=[A,zeros(wXh_Dry)]
B=ones(size(A)); C=[A;B]  ;
D=repmat(C,3,1);
D=[B;B;D;B];
imshow(D,[]) ; % monitor the pore space
Channel2D=D;
Len_Channel_2D=size(Channel2D,1); % Length
Width=size(Channel2D,2); % should not be hod
Channel_2D_half_Width=Width/2,

% test without obstacles (i.e. 2D channel & no obstacles)

if (Pois_test)
%over-writes the definition of the pore space
clear Channel2D
Len_Channel_2D=18,
Channel_2D_half_Width=4; Width=Channel_2D_half_Width*2;
Channel2D=ones(Len_Channel_2D,Width); % define wet area
%Channel2D(6:12,6:8)=0; % put fluid obstacle
imshow(Channel2D,[]);
end

[Nr Mc]=size(Channel2D); % Number rows and Munber columns

% FLUID PROPERTIES
% physical properties
cs2=1/3; %
cP_visco=0.5; % [cP] 1 CP Dinamic water viscosity 20 C
density=1.; % fluid density
Lky_visco=cP_visco/density; % lattice kinematic viscosity
omega=(Lky_visco/cs2+0.5).^-1; %  omega: relaxation frequency
%Lky_visco=cs2*(1/omega - 0.5) , % lattice kinematic viscosity
%dPdL= Pressure / dL;% External pressure gradient [atm/cm]

uy_fin_max=-0.2;
%dPdL = abs( 2*Lky_visco*uy_fin_max/(Channel_2D_half_Width.^2) );
dPdL=-0.0125;
uy_fin_max=dPdL*(Channel_2D_half_Width.^2)/(2*Lky_visco); % Poiseuille Gradient;
% max poiseuille final  velocity on the flow profile
uy0=-0.001; ux0=0.0001; %  linear vel .. inizialization

%
% uy_fin_max=-0.2; % max poiseuille final  velocity on the flow profile
% omega=0.5, cs2=1/3; % omega: relaxation frequency
% Lky_visco=cs2*(1/omega - 0.5) , % lattice kinematic viscosity
% dPdL = abs( 2*Lky_visco*uy_fin_max/(Channel_2D_half_Width.^2) ); % Poiseuille Gradient;
%

uyf_av=uy_fin_max*(2/3);; % average fluid velocity on the profile

x_profile=([-Channel_2D_half_Width:+Channel_2D_half_Width-1]+0.5);
uy_analy_profile=uy_fin_max.*(1-  ( x_profile /Channel_2D_half_Width).^2 ); % analytical velocity profile

av_vel_t=1.e+10; % inizialization (t=0)
%PixelSize= 5; % [Microns]
%dL=(Nr*PixelSize*1.0E-4); % sample hight [cm]

%
% EXPERIMENTAL SET-UP
% inlet and outlet buffers
inb=2, oub=2; % inlet and outlet buffers thickness
% add fluid at the inlet (top) and outlet (down)
inlet=ones(inb,Mc); outlet=ones(oub,Mc);
Channel2D=[ [inlet]; Channel2D ;[outlet] ] ; % add flux in and down (E to W)
[Nr Mc]=size(Channel2D); % update size
% boundaries related to the experimental set up
wb=2; % wall thickness
Channel2D=[zeros(Nr,wb), Channel2D , zeros(Nr,wb)]; % add walls (no fluid leak)
[Nr Mc]=size(Channel2D); % update size
uy_analy_profile=[zeros(1,wb), uy_analy_profile, zeros(1,wb) ] ; % take into account walls
x_pro_fig=[[x_profile(1)-[wb:-1:1]], [x_profile, [1:wb]+x_profile(end)] ];

% Figure plots analytical parabolic profile
figure(20), plot(x_pro_fig,uy_analy_profile,'-'), grid on,
title('Analytical parab. profile for Poiseuille planar flow in a channel')

% VISUALIZE PORE SPACE & FLUID OSTACLES & MEDIAL AXIS
figure, imshow(Channel2D); title('Vassel geometry');
Channel2D=logical(Channel2D);
% obstacles for Bounce Back ( in front of the grain)
Obstacles=bwperim(Channel2D,8); % perimeter of the grains for bounce back Bound.Cond.
Obstacles([1:inb,Nr-oub:Nr],[wb+2:Mc-wb-1])=0;

figure, imshow(Obstacles); title(' Fluid obstacles (in the fluid)' );
%
Medial_axis=bwmorph(Channel2D,'thin',Inf); %
figure, imshow(Medial_axis); title('Medial axis');
%figure(10) % used to visualize evolution of rho
%figure(11) % used to visualize ux
%figure(12) % used to visualize uy (i.e. top -> down)

% porosity
porosity=nnz(Channel2D==1)/(Nr*Mc)% porosity
porosity=nnz(Channel2D( :,[(wb+1):(Mc-wb)] )==1)./(Nr*(Mc-2*wb))
% exclude the wall from counting

% INDICES
% Wet locations etc.
[iabw1 jabw1]=find(Channel2D==1); % indices i,j, of active lattice locations i.e. pore
lena=length(iabw1); % number of active location i.e. of pore space lattice cells
ija= (jabw1-1)*Nr+iabw1; % equivalent single index (i,j)->> ija for active locations
% absolute (single index) position of the obstacles in for bounce back in Channel2D
% Obstacles
[iobs jobs]=find(Obstacles);lenobs=length(iobs); ijobs= (jobs-1)*Nr+iobs; % as above
% Medial axis of the pore space
[ima jma]=find(Medial_axis); lenma=length(ima);  ijma= (jma-1)*Nr+ima; % as above
% Internal wet locations : wet & ~obstables
% (i.e. internal wet lattice location non in contact with dray locations)
[iawint jawint]=find(( Channel2D==1 & ~Obstacles)); % indices i,j, of active lattice locations
lenwint=length(iawint); % number of internal (i.e. not border) wet locations
ijaint= (jawint-1)*Nr+iawint; % equivalent singl
NxM=Nr*Mc;

% Look-up table for direct acess to the position of the ija in a wet only
% array , i.e. from a i,j position in an array (full) to the index in an arry of wet
% positions only, eg the position 5 is the 3 because 2 pixels are dry

%ija =[2 3 5 9 13]; % wet positions
Lut(:,1)=[1:1:ija(end)]'; Lut(ija,2)=[1:1:length(ija)];
% Lut use ijac=9
% from current ija (dry & wet ) to its position in an array of only wet
% Lut(ijac,2) ... % alternative find(ija==ijac)

% DIRECTIONS: E N W S NE NW SW SE ZERO (ZERO:Rest Particle)
%    y^
%  6 2 5           ^         NW  N  NE
%  3 9 1 ... +x-> +y         W   RP  E
%  7 4 8                     SW  S  SE
%   -y
% x & y components of velocities , +x is to est , +y is to nord
East=1; North=2; West=3; South=4; NE=5; NW=6; SW=7; SE=8; RP=9;
N_c=9 ; % number of directions
% versors D2Q9
C_x=[1 0 -1  0 1 -1 -1  1 0];
C_y=[0 1  0 -1 1  1 -1 -1 0]; C=[C_x;C_y]

% BOUNCE BACK SCHEME
% after collision the fluid elements densities f are sent back to the
% lattice node they come from with opposite direction
% indices opposite to 1:8 for fast inversion after bounce
ic_op = [3 4 1 2 7 8 5 6]; %   i.e. 4 is opposite to 2 etc.

% PERIODIC BOUNDARY CONDITIONS - reinjection rules
yi2=[Nr , 1:Nr , 1]; % this definition allows implemening Period Bound Cond
%yi2=[1, Nr , 2:Nr-1 , 1,Nr]; % re-inj the second last to as first
% directional weights (density weights)
w0=16/36. ; w1=4/36. ; w2=1/36.;
W=[ w1 w1 w1 w1 w2 w2 w2 w2 w0];
%c constants (sound speed related)
cs2=1/3; cs2x2=2*cs2; cs4x2=2*cs2.^2;
f1=1/cs2; f2=1/cs2x2; f3=1/cs4x2;
f1=3.; f2=4.5; f3=1.5; % coef. of the f equil.

% declarative statemets
f=zeros(lena,N_c); % array of fluid density distribution
feq=zeros(lena,N_c); % f at equilibrium
rho=ones(lena,1); % macro-scopic density
temp1=zeros(lena,1);
ux=zeros(lena,1);   uy=zeros(lena,1); uyout=zeros(lena,1);  % dimensionless velocities
uxsq=zeros(lena,1); uysq=zeros(lena,1);   usq=zeros(lena,1);  % higher degree velocities
% for normal 2d visualization
uy1=zeros(Nr,Mc);  ux1=zeros(Nr,Mc);

% initialization arrays : start values in the wet area
for ia=1:lena % stat values in the active cells only ; 0 outside
    f(ia,:)=1/9; % uniform density distribution for a start
end
uy(:)=uy0; ux(:)=ux0; % initialize fluid velocities
rho(:)=density;

% EXTERNAL (Body) FORCES e.g. inlet pressure or inlet-outlet gradient
% directions: E N W S NE NW SW SE ZERO
force = -dPdL*(1/6)*1*[0 -1 0 1 -1 -1 1  1  0]'; %;
%...                   E  N E S NE NW SW SE RP ...
% the pressure pushes the fluid down i.e. N to S

% While .. MAIN TIME EVOLUTION LOOP
StopFlag=false; % i.e. logical(0)
Max_Iter=3000; % max allowed number of iteration
Check_Iter=1; Output_Every=8; % frequency of check & output
Cur_Iter=0; % current iteration counter inizialization
toler=1.0e-8; % tollerance to declare convegence
Cond_path=[]; % recording values of the convergence criterium
density_path=[]; % recording aver. density values for convergence
end % ends if restart

if(Restart==true)
 StopFlag=false;  Max_Iter=Max_Iter+3000; toler=1.0e-12;
end

while(~StopFlag)
    Cur_Iter=Cur_Iter+1 % iteration counter update

    % density and moments
    rho=sum(f,2); % density

    if Cur_Iter >1 % use inizialization ux uy to start
        % Moments ... Note:C_x(9)=C_y(9)=0
       % ux=zeros(lena,1); uy=zeros(lena,1);
        %for ic=1:N_c-1;
         %   ux = ux + C_x(ic).*f(:,ic) ; uy = uy + C_y(ic).*f(:,ic)  ;
       % end

        uy=f(:,2) +f(:,5)+f(:,6)-f(:,4)-f(:,7)-f(:,8); % in short !
        ux=f(:,1) +f(:,5)+f(:,8)-f(:,3)-f(:,6)-f(:,7); % in short !

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ux(:)=ux(:)./rho(:); uy(:)=uy(:)./rho(:);
    uxsq(:)=ux(:).^2; uysq(:)=uy(:).^2; usq(:)=uxsq(:)+uysq(:); %

    % weighted densities : rest particle, principal axis, diagonals
    rt0 = w0.*rho; rt1 = w1.*rho; rt2 = w2.*rho;

    % Equilibrium distribution
    % main  directions ( + cross)
    feq(:,1)= rt1(:) .*(1 +f1*ux(:) +f2*uxsq(:) -f3*usq(:));
    feq(:,2)= rt1(:) .*(1 +f1*uy(:) +f2*uysq(:) -f3*usq(:));
    feq(:,3)= rt1(:) .*(1 -f1*ux(:) +f2*uxsq(:) -f3*usq(:));
    %feq(ija+NxM*(3)=f(ija)-2*rt1(ija)*f1.*ux(ija); % much faster... !!
    feq(:,4)= rt1(:) .*(1 -f1*uy(:) +f2*uysq(:) -f3*usq(:));

    % diagonals (X diagonals) (ic-1)
    feq(:,5)= rt2(:) .*(1 +f1*(+ux(:)+uy(:)) +f2*(+ux(:)+uy(:)).^2 -f3.*usq(:));
    feq(:,6)= rt2(:) .*(1 +f1*(-ux(:)+uy(:)) +f2*(-ux(:)+uy(:)).^2 -f3.*usq(:));
    feq(:,7)= rt2(:) .*(1 +f1*(-ux(:)-uy(:)) +f2*(-ux(:)-uy(:)).^2 -f3.*usq(:));
    feq(:,8)= rt2(:) .*(1 +f1*(+ux(:)-uy(:)) +f2*(+ux(:)-uy(:)).^2 -f3.*usq(:));
    % rest particle (.) ic=9
    feq(:,9)= rt0(:) .*(1 - f3*usq(:));

    %Collision (between fluid elements)omega=relaxation frequency
    f=(1.-omega).*f + omega.*feq;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %add external body force due to the pressure gradient prop. to dPdL
    for ic=1:N_c;%-1
       % for ia=1:lena
           % i=iabw1(ia);  j=jabw1(ia);
            % if Obstacles(i,j)==0 % the i,j is not aderent to the boundaries
            f(:,ic)= f(:,ic) + force(ic);
            % end
       % end
    end

    % % STREAM
    % Forward Propagation step & % Bounce Back (collision fluid with obstacles)
    %f(:,9) = f(:,9); % Rest element do not move

    feq = f; % temp storage of f in feq
        for ic=1:1:N_c-1, % select velocity layer
        C_yic = C_y(ic); C_xic = C_x(ic);
        ic2=ic_op(ic); % selects the layer of the velocity opposite to ic for BB
        temp1=feq(:,ic); %

        % from wet location that are NOT on the border to other wet locations
        for ia=1:1:lenwint % number of internal (i.e. not border) wet locations
            i=iawint(ia);  j=jawint(ia);  % so that we care for the wet space only !
            i2 = i+C_yic; j2 = j+C_xic; % Expected final locations to move
            i2=yi2(i2+1); % i2 corrected for PBC when necessary (flow out re-fed to inlet)
            % i.e the new position (i2,j2)==ia2 is sure another wet location
            % therefore normal propagation from (i,j)==ijaint to (i2,j2)==ia2 on layer ic
            ia2=(j2-1)*Nr+i2; iarec2=Lut(ia2,2) ; % use lut (look up table)
            iarec1=Lut(ijaint(ia),2) ;
            f(iarec2,ic)=temp1(iarec1); % copy
        end ; % i and j single loop

        % from wet locations that ARE on the border of obstacles
        for ia=1:1:lenobs % wet border locations
            i=iobs(ia);  j=jobs(ia);  % so that we care for the wet space only !
            i2 = i+C_yic; j2 = j+C_xic; % Expected final locations to move
            i2=yi2(i2+1); % i2 corrected for PBC
            iarec1=Lut(ijobs(ia),2) ;
            if( Channel2D(i2,j2) ==0 ) % i.e the new position (i2,j2) is dry
               f(iarec1,ic2) =temp1(iarec1); % invert direction: bounce-back in the opposite direction ic2
            else % otherwise, normal propagation from (i,j) to (i2,j2) on layer ic
               ia2=(j2-1)*Nr+i2; iarec2=Lut(ia2,2) ;
               f(iarec2,ic)=temp1(iarec1); %
            end ; % b.b. and propagations

        end ; % i and j single loop
        % special treatment for Corners
        %   f(1,wb+1,ic)=temp1(Nr,Mc-wb);      f(1,Mc-wb,ic)=temp1(Nr,wb+1);
        %   f(Nr,wb+1,ic)=temp1(1,Mc-wb);      f(Nr,Mc-wb,ic)=temp1(1,wb+1);

    end ; %  for ic direction

    % ends of Forward Propagation step &  Bounce Back Sections

    % re-calculate  uy as uyout for convergence
    rho=sum(f,2); % density
    % check velocity
%     uyout= zeros(lena,1);
%     for ic=1:N_c-1;
%         uyout= uyout + C_y(ic).*f(:,ic) ; % flow dim.less velocity out
%     end
%    % uyout(ija)=uyout(ija)./rho(ija); % from momentum to velocity

    % Convergence check on velocity values
    if (mod(Cur_Iter,Check_Iter)==0) ; % check for convergence every 'Check_Iter' iterations

        % variables monitored
        % mean density and
        vect=rho(:); vect=vect(:);
        cur_density=mean(vect);
        % mean 'interstitial' velocity
        % uy(ija)=uy(ija)/rho(ija); ?
        vect=uy(:); av_vel_int= mean(vect)  ; % seepage velocity (in the wet area)
        % on the whole cross-sectional area of flow (wet + dry)
        av_vel_int=av_vel_int*porosity, % av. vel. on the wet + dry area
        %av_vel_int=mean2(uy),
        av_vel_tp1 = av_vel_int;
        Condition=abs( abs(av_vel_t/av_vel_tp1 )-1), % should --> 0

        Cond_path=[Cond_path, Condition]; % records the convergence path (value)
        density_path=[density_path, cur_density];
        %
        av_vel_t=av_vel_tp1; % time t & t+1

        if (Condition < toler) | (Cur_Iter > Max_Iter)
            StopFlag=true;
            display( 'Stop iteration: Convergence met or iteration exeeding the max allowed' )
            display( ['Current iteration: ',num2str(Cur_Iter),...
                ' Max Number of iter: ',num2str(Max_Iter)] )
            break % Terminate execution of WHILE .. exit the time evolution loop.

        end    % if(Condition < toler

    end

    if (mod(Cur_Iter,Output_Every)==0) ;  % Output from loop every ...
        %if (Cur_Iter>60) ;  % Output from loop every ...

       % rho=sum(f,2); % density
       % figure(10); imshow(rho,[0.1 0.9]); title(' rho'); % visualize density evolution
        uy1(ija)=uy;  ux1(ija)=ux;
        %figure(11); imshow(ux1,[ ]); title(' ux' ); % visualize fluid velocity horizontal
        %figure(12); imshow(-uy1,[ ]); title(' uy' ); % visualize fluid velocity down
        %figure(14), imshow(-uyout,[]), title('uyout'); % vis vel flow out
        up=2; % linear section to visualize up from the lower row
        figure(15), hold off, feather(ux1(Nr-up,:),uy1(Nr-up,:)),
        figure(15), hold on , plot(uy_analy_profile,'r-')
        title('Analytical (i.e. no obstacles) vs LB calculated (Blue), fluid velocity profile')
        pause(1); % time given to visualize properly

    end % every

   % pause(1);

end %  End main time Evolution Loop

% Output & Draw after the end of the time evolution

figure, plot(Cond_path(2:end)); title('convergence path')
%figure, plot(density_path(2:end)); title('density convergence path')
figure, plot( [uy1(Nr-up,:)-uy_analy_profile] ); title('difference : LB - Analytical solution')

toc

% Permeability K

K_Darcy_Porous_Sys= (av_vel_int*porosity)/dPdL*Lky_visco ,

K_Analy_2D_Channel=(Width^2)/12

repo: Aero-matlab


text:
function e = end(v,i,n)
% overloaded end function

if n==1
  e = numel(v.x);
else
  e = size(v,i);
end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Simulation_close_tool:
%   Set the tool_open variable to zero. When drawn using the drawrobot3d function
%   the tool will be drawn with the coordinates defined as open.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function simulation_close_tool

global configuration robot

robot.tool.tool_open=0;

drawrobot3d(robot, robot.q);

repo: ARTE


text:
function weight = read_filter(file_name);
FID = fopen(file_name,'r');
tmp = textscan(FID,'%f ',2);
filter_len = tmp{1}(1);
N_ch = tmp{1}(2);
for i=1:N_ch,
    tmp = textscan(FID,'%f ', filter_len);
    weight(:,i) = tmp{1};
end
fclose(FID);

author: singaxiong

repo: SignalGraph


text:
% Fig. 6.30   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

num=[1 1];
den=conv([1 0],[0.1 -1]);
w=logspace(-1,2,100);
[m,p]=bode(num,den,w);
subplot(2,1,1),loglog(w,m,'-',w,ones(size(w)),'-');
grid;
ylabel('Magnitude');
title('Fig. 6.30 Bode Plot for G=(s+1)/[s(s/10 -1)] (a) magnitude');
subplot(2,1,2)
if p>0, p=p-360; end      % Matlab quadrant control varies in different versions
semilogx(w,p,'-',w,-180*ones(size(w)),'-');
axis([.1 100 -270 -90])
grid;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('(b) phase');

author: Sable

repo: mcbench-benchmarks


text:
function [ colptr, rowind ] = hb_structure_read ( input_unit, ncol, mxtype, ...
  nnzero, neltvl, ptrcrd, ptrfmt, indcrd, indfmt )

%*****************************************************************************80
%
%% HB_STRUCTURE_READ reads the structure of an HB matrix.
%
%  Discussion:
%
%    The user should already have opened the file, and positioned it
%    to just after the header records.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the unit from which data is read.
%
%    Input, integer NCOL, the number of columns.
%
%    Input, character ( len = 3 ) MXTYPE, the matrix type.
%    First character is R for Real, C for complex, P for pattern only.
%    Second character is S for symmetric, U for unsymmetric, H for
%      Hermitian, Z for skew symmetric, R for rectangular.
%    Third character is A for assembled and E for unassembled
%      finite element matrices.
%
%    Input, integer NNZERO.  In the case of assembled sparse matrices,
%    this is the number of nonzeroes.  In the case of unassembled finite
%    element matrices, in which the right hand side vectors are also
%    stored as unassembled finite element vectors, this is the total
%    number of entries in a single unassembled right hand side vector.
%
%    Input, integer NELTVL, the number of finite element matrix entries,
%    set to 0 in the case of assembled matrices.
%
%    Input, integer PTRCRD, the number of pointer records.
%
%    Input, character ( len = 16 ) PTRFMT, the format for reading pointers.
%
%    Input, integer INDCRD, the number of index records.
%
%    Input, character ( len = 16 ) INDFMT, the format for reading indices.
%
%    Output, integer COLPTR(NCOL+1), COLPTR(I) points to the location of
%    the first entry of column I in the sparse matrix structure.
%
%    Output, integer ROWIND(NNZERO) or ROWIND(NELTVL), the row index of
%    each item.
%
  [ p, code, w, m ] = s_to_format ( ptrfmt );

  if ( mxtype(3) == 'A' )
    line_num = 1 + floor ( ( ( ncol + 1 ) - 1 ) / p );
  else
    line_num = 1 + floor ( ( ( ncol     ) - 1 ) / p );
  end

  jhi = 0;

  for i = 1 : line_num
    line = fgetl ( input_unit );
    jlo = jhi + 1;
    if ( mxtype(3) == 'A' )
      jhi = min ( jlo + p - 1, ncol + 1 );
    else
      jhi = min (  jlo + p - 1, ncol );
    end

    khi = 0;

    for j = jlo : jhi
      klo = khi + 1;
      khi = min ( klo + w - 1, length ( line ) );
      s = line(klo:khi);
      colptr(j) = s_to_i4 ( s );
    end

  end

  if ( mxtype(3) == 'A' )

    [ p, code, w, m ] = s_to_format ( indfmt );

    line_num = 1 + floor ( ( nnzero - 1 ) / p );

    jhi = 0;

    for i = 1 : line_num
      line = fgetl ( input_unit );
      jlo = jhi + 1;
      jhi = min ( jlo + p - 1, nnzero );

      khi = 0;

      for j = jlo : jhi
        klo = khi + 1;
        khi = min ( klo + w - 1, length ( line ) );
        s = line(klo:khi);
        rowind(j) = s_to_i4 ( s );
      end

    end

  elseif ( mxtype(3) == 'E' )

    [ p, code, w, m ] = s_to_format ( indfmt );
    number = colptr(ncol) - colptr(1);
    line_num = 1 + floor ( ( number - 1 ) / p );

    jhi = 0;

    for i = 1 : line_num
      line = fgetl ( input_unit );
      jlo = jhi + 1;
      jhi = min ( jlo + p - 1, number );

      khi = 0;

      for j = jlo : jhi
        klo = khi + 1;
        khi = min ( klo + w - 1, length ( line ) );
        s = line(klo:khi);
        rowind(j) = s_to_i4 ( s );
      end

    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'HB_STRUCTURE_READ - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of MXTYPE(3).\n' );
    return

  end

  return
end


text:
function F = conj(F)
%CONJ   Complex conjugate of a CHEBFUN3V.
%   CONJ(F) returns the complex conjugate of F. For a complex F, we have
%   CONJ(F) = REAL(F) - i*IMAG(F).
%
% See also CHEBFUN3V/REAL and CHEBFUN3V/IMAG.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty(F) )
    return
end

% Take conj part of each component:
F.components = cellfun(@conj, F.components, 'UniformOutput', false);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
classdef CAEAD < ALGORITHM
% <multi> <real/integer/label/binary/permutation> <constrained>
% Dual-population evolutionary algorithm based on alternative evolution and degeneration
% type --- 1 --- Type of operator (1. DE 2. GA)

%------------------------------- Reference --------------------------------
% J. Zou, R. Sun, S. Yang, and J. Zheng, A dual-population algorithm based
% on alternative evolution and degeneration for solving constrained multi-
% objective optimization problems, Informaction Scinece, 2021, 239: 89-102.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB Platform
% for Evolutionary Multi-Objective Optimization [Educational Forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

   methods
       function main(Algorithm,Problem)
           type = Algorithm.ParameterSet(1);

           %% Generate random population
           Population1 = Problem.Initialization();
           Population2 = Problem.Initialization();
           Fitness1    = CalFitness(Population1.objs,Population1.cons,0);
           Fitness2    = CalFitness(Population2.objs,Population2.cons,1e6);

           min_epsilon      = 1e-4;
           change_threshold = 1e-2;
           max_change       = 1;
           epsilon_k        = 1e8;
           tao              = 0.05;
           max_ep           = 0;
           gen              = 1;
           stage            = false;

           %% Optimization
           while Algorithm.NotTerminated(Population1)
               pop_cons2      = Population2.cons;
               cv2            = overall_cv(pop_cons2);
               population     = [Population2.decs,Population2.objs,cv2];
               Objvalues(gen) = sum(sum(Population2.objs,1));
               ep(gen)        =  epsilon_k;
               if type == 1
                   MatingPool1 = TournamentSelection(2,2*Problem.N,Fitness1);
                   MatingPool2 = TournamentSelection(2,2*Problem.N,Fitness2);
                   Offspring1  = OperatorDE(Problem,Population1,Population1(MatingPool1(1:end/2)),Population1(MatingPool1(end/2+1:end)));
                   Offspring2  = OperatorDE(Problem,Population2,Population2(MatingPool2(1:end/2)),Population2(MatingPool2(end/2+1:end)));
               elseif type == 2
                   MatingPool1 = TournamentSelection(2,Problem.N,Fitness1);
                   MatingPool2 = TournamentSelection(2,Problem.N,Fitness2);
                   Offspring1  = OperatorGAhalf(Problem,Population1(MatingPool1));
                   Offspring2  = OperatorGAhalf(Problem,Population2(MatingPool2));
               end
               [FrontNo2,~] = NDSort(Population2.objs,size(Population2.objs,1));
               NC2 = size(find(FrontNo2==1),2);
               if gen ~= 1
                   max_change = abs(Objvalues(gen)-Objvalues(gen-1));
               end
               if max_change <= change_threshold &&NC2 == Problem.N && stage == false
                   epsilon_k = max(population(:,end),[],1);
                   stage = true;
               end
               Offspring3 = [];
               if stage == true
                   if type == 1
                       Offspring3 = OperatorDE(Problem,Population1,Population2(MatingPool2(1:end/2)),Population2(MatingPool2(end/2+1:end)));
                   elseif type == 2
                       for i=1:Problem.N/2
                           Offtemp = OperatorGAhalf(Problem,[Population1(MatingPool1(i)),Population2(MatingPool2(i))]);
                           Offspring3 = [Offspring3,Offtemp];
                       end
                   end
               end
               if stage == true
                   [stage,epsilon_k] =  update_epsilon(stage,tao,epsilon_k,max_ep,min_epsilon);
               end
               if epsilon_k < 9e5
                   max_ep = max(max_ep,epsilon_k);
               end
               [Population1,Fitness1] = EnvironmentalSelection([Population1,Offspring1,Offspring2,Offspring3],Problem.N,true,0);
               [Population2,Fitness2] = EnvironmentalSelection([Population2,Offspring2],Problem.N,false,epsilon_k);
               gen = gen+1;
           end
       end
   end
end

function result = overall_cv(cv)
    cv(cv <= 0) = 0;cv = abs(cv);
    result = sum(cv,2);
end

function [stage,result] = update_epsilon(stage,tao,epsilon_k,epsilon_0,min_epsilon)
    if epsilon_k > min_epsilon
        result = (1 - tao) * epsilon_k;
        stage = true;
    else
        result = epsilon_0;
        stage = false;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function hit_node(varargin)
%Function for generating Another metabolite-metabolite interaction network
%if the node is clicked in the figure that produced by createMetIntrcNetwork.
%Right click and left click has different properties.
%Left click :Generate sub-metabolite-metabolite network from the created figure,
%            this functionality was added for better looking at the created network
%            and showing flux values on edges lines.
%
%Right click :Generate metabolite metabolite network from model, this
%             functionality were added for creating metabolite-metabolite network using
%             clicked metabolite and model.
%The produced figure from main figure also has same property with main
%figure, and it is clickable too.

    f = varargin{1}.Parent.Parent;
    if (strcmp(f.SelectionType, 'normal') | strcmp(f.SelectionType, 'alt'))
        metofInt=varargin{1}.String;
        adjncMtrx=varargin{3};
        mets=varargin{4};
        model=varargin{5} ;
        excludedMets=varargin{7};
        fluxes=varargin{6};
        Graphtitle=varargin{8};
        nodecolour=varargin{9};
        Hnodecolour=varargin{10};
        fcont=varargin{11};
        scaleMin=varargin{12};
        scaleMax=varargin{13};
        nodeSize=varargin{14};
        HnodeSize=varargin{15};
        arrowSize=varargin{16};
        threshold=varargin{17};
        excNodesWithDeg=varargin{18};

% Define function of left click
        if strcmp(f.SelectionType, 'normal')
            idx=ismember(mets,metofInt);
            newadjMtrx=zeros(size(adjncMtrx));
            newadjMtrx(idx,:)=adjncMtrx(idx,:);
            newadjMtrx(:,idx)=adjncMtrx(:,idx);
            rows=any(newadjMtrx,1);
            cols=any(newadjMtrx,2);
            indofMets=rows'|cols;
            metnames=mets(indofMets);
            newadjMtrx=newadjMtrx(indofMets,indofMets);
            A=newadjMtrx;
            G = digraph(A,metnames);
            ttl=['Sub-', ' ',Graphtitle];
        end
% Define function of right click
        if strcmp(f.SelectionType, 'alt')
            IndexesofMets=find(ismember(model.mets,metofInt));
            metMatrix=~ismember(model.S(IndexesofMets,:),0);
            Rxns=model.rxns(any(metMatrix,1));
            FluxRes=fluxes(ismember(model.rxns,Rxns));
            metmatrix=~ismember(model.S(:,any(metMatrix,1)),0);
            metnames=model.mets(any(metmatrix,2));
            adjunc_mat=model.S(any(metmatrix,2),any(metMatrix,1));
            FluxMatrix  =FluxRes + 1e-6;
            FluxMatrix = repmat(FluxMatrix',size(adjunc_mat,1),1);
            adjunc_mat=adjunc_mat.*FluxMatrix;
            leftMatrix=adjunc_mat;
            leftMatrix(leftMatrix > 0)=0;
            rightMatrix=adjunc_mat;
            rightMatrix(rightMatrix < 0)=0;
            rightMatrix(rightMatrix > 0)=1;
            metMatrix=leftMatrix*rightMatrix';
            A=metMatrix*(-1);
            idx=ismember(metnames,metofInt);
            newadjMtrx=zeros(size(A));
            newadjMtrx(idx,:)=A(idx,:);
            newadjMtrx(:,idx)=A(:,idx);
            rows=any(newadjMtrx,1);
            cols=any(newadjMtrx,2);
            indofMets=rows'|cols;
            metnames=metnames(indofMets);
            newadjMtrx=newadjMtrx(indofMets,indofMets);
            A=newadjMtrx;
            G = digraph(A,metnames);
            ttl=Graphtitle;
        end

        G = rmnode(G,excludedMets);
        G.Edges.Weight(G.Edges.Weight < 1e-4)=1e-6;
        G.Edges.LWidths = (G.Edges.Weight-min(G.Edges.Weight))/(max(G.Edges.Weight)-min(G.Edges.Weight))+0.00001;
        G.Edges.LWidths((isnan(G.Edges.LWidths)))=1e-6;

        if threshold~=1e-7
            edesBelowTresholdIdx=find(G.Edges.Weight<threshold);
            G=rmedge(G,edesBelowTresholdIdx);
        end

        nodesIndegree=indegree(G);
        nodesoutdegree=outdegree(G);
        totalDegree=nodesIndegree+nodesoutdegree;
        G=rmnode(G,find(totalDegree<1));

        figure;
        hold on
        h=plot(G,'MarkerSize',(nodeSize+5),'NodeColor',nodecolour,'ArrowSize',(arrowSize+5),'NodeLabelMode','auto');
        layout(h,'layered','Direction','right')
        highlight(h,metofInt,'NodeColor',Hnodecolour,'MarkerSize',(HnodeSize+5));
        nl = h.NodeLabel;
        h.NodeLabel = '';
        xd = get(h, 'XData');
        yd = get(h, 'YData');
        title([metofInt,' ', 'Centred', ' ',ttl],'Interpreter', 'none');
        txt=text(xd, yd, nl, 'FontSize',8, 'FontWeight','bold', 'HorizontalAlignment','center', 'VerticalAlignment','middle');
        set(txt,'Interpreter', 'none');
        set(txt,'ButtonDownFcn',{@hit_node,A,metnames,model,fluxes,excludedMets,Graphtitle,nodecolour,Hnodecolour,fcont,scaleMin,scaleMax,nodeSize,HnodeSize,arrowSize,threshold,excNodesWithDeg});
        if  fcont==1
            h.LineWidth=1;
        elseif ~any(isnan(G.Edges.LWidths))
            h.EdgeLabel=G.Edges.Weight;
            h.LineWidth=G.Edges.LWidths*2.5;
            colormap jet(10)
            h.EdgeCData=G.Edges.Weight;
            hcb=colorbar;
            colorTitleHandle = get(hcb,'Title');

        else
            h.LineWidth=1;
            colormap jet(10)
            h.EdgeCData=G.Edges.Weight;
            hcb=colorbar;
            colorTitleHandle = get(hcb,'Title');
        end

        if  fcont==0
        if ~(scaleMax==1e-6 && scaleMin==0)
            caxis([scaleMin scaleMax])
            titleBar = ['Scaled Fluxes',' ', '(', num2str(scaleMin),' - ',num2str(scaleMax),')'];
            set(colorTitleHandle ,'String',titleBar,'FontWeight','Bold');
        else
            set(colorTitleHandle ,'String','Fluxes','FontWeight','Bold');
        end
        end

        set(gca,'XTickLabel',{' '});
        set(gca,'YTickLabel',{' '});
        set(gca,'YTick',[]);
        set(gca,'XTick',[]);
        set(gca,'XColor', 'none','YColor','none');
    end

end

author: opencobra

repo: cobratoolbox


text:
function triangle_ncc_rule_test04 ( )

%*****************************************************************************80
%
%% TEST04 tests TRIANGLE_NCC_RULE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST04\n' );
  fprintf ( 1, '  TRIANGLE_NCC_RULE returns the points and weights of\n' );
  fprintf ( 1, '  an NCC rule for the unit triangle.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This routine uses those rules to estimate the\n' );
  fprintf ( 1, '  integral of monomomials in the unit triangle.\n' );

  rule_num = triangle_ncc_rule_num ( );

  area = 0.5;

  for a = 0 : 10

    for b = 0 : 10 - a
%
%  Multiplying X^A * Y^B by COEF will give us an integrand
%  whose integral is exactly 1.  This makes the error calculations easy.
%
      coef = ( a + b + 2 ) * ( a + b + 1 );
      for i = 1 : b
        coef = coef * ( a + i ) / i;
      end

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Integrate %f * X^%d * Y^%d\n', coef, a, b );
      fprintf ( 1, '\n' );
      fprintf ( 1, '      Rule       QUAD           ERROR\n' );
      fprintf ( 1, '\n' );

      for rule = 1 : rule_num

        order_num = triangle_ncc_order_num ( rule );

        [ xy, w ] = triangle_ncc_rule ( rule, order_num );

        quad = 0.0;

        for order = 1 : order_num

          x = xy(1,order);
          y = xy(2,order);

          if ( a == 0 & b == 0 )
            value = coef;
          elseif ( a == 0 & b ~= 0 )
            value = coef * y^b;
          elseif ( a ~= 0 & b == 0 )
            value = coef * x^a;
          elseif ( a ~= 0 & b ~= 0 )
            value = coef * x^a * y^b;
          end

          quad = quad + w(order) * value;

        end

        quad = area * quad;

        exact = 1.0;
        err = abs ( exact - quad );

        fprintf ( 1, '  %8d  %14f  %14f\n', rule, quad, err );

      end

    end

  end

  return
end


text:
function c = cellRemoveEmpty(inCell)
%Remove empty entries from a cell array
%
%   c = cellRemoveEmpty(inCell)
%
% GB/BW
% See removeListElement for another way to do this.
%

if isempty(inCell), c = {}; return; end
n = length(inCell);
emptyList = [];
for ii=1:n
    if isempty(inCell{ii}), emptyList = [emptyList,ii]; end
end
if isempty(emptyList), c = inCell; return; end
c = {};
for ii=setdiff(1:n,emptyList)
    c{end+1} = inCell{ii};
end
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef TestEllipse2Poly
    %TestEllipse2Poly

    methods (Static)
        function test_1
            pts = cv.ellipse2Poly([64,64], [20,10]);
            validateattributes(pts, {'numeric'}, ...
                {'2d', 'size',[NaN 2], 'integer'});
        end

        function test_2
            pts = cv.ellipse2Poly([64,64], [20,10], 'Angle',30, ...
                'StartAngle',15, 'EndAngle',200, 'Delta',2);
            validateattributes(pts, {'numeric'}, ...
                {'2d', 'size',[NaN 2], 'integer'});
        end

        function test_3
            pts = cv.ellipse2Poly([64,64], [20,10], 'DoublePrecision',true);
            validateattributes(pts, {'numeric'}, ...
                {'2d', 'size',[NaN 2], 'real'});
        end

        function test_error_argnum
            try
                cv.ellipse2Poly();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [message,result,solverTime] = runcontest(flagVisualize,charFlag,whichBoards)
%RUNCONTEST Test an entry.
%   [MESSAGE,RESULTS,TIME] = RUNCONTEST runs the file solver.m against all
%   the problems defined in testsuite_sample.mat. MESSAGE returns a summary
%   of the testing, RESULTS measures how well the entry solved the problem,
%   and TIME measures the time the entry took to compute the answer.
%
%   RUNCONTEST(true) graphically visualize the results.
%
%   RUNCONTEST(true/false, true) displays numbers as their ASCII
%   equivalents
%
%   RUNCONTEST(true/false, true/false, DOBOARDS) runs the file solver.m
%   against the problems enumerated in the vector DOBOARDS as defined in
%   the testsuite.

% Copyright 2011 The MathWorks, Inc.

% The MATLAB Contest Team
% April 2011

testSuiteFile = 'testsuite_sample.mat';
tests = load(testSuiteFile,'testsuite');

% Argument parsing.
if nargin < 3
    whichBoards = 1:numel(tests.testsuite); % do all boards in test suite
end
if nargin < 2
    charFlag = false;  % display numbers as numbers
end
if nargin < 1
    flagVisualize = false;  % no visualization
end

numBoards = numel(whichBoards);
solverTime = zeros(numBoards,1);
solverScore = zeros(numBoards,1);

if flagVisualize
    fh = figure('Color', [1 1 1]);
end

for n = 1:numBoards
    [words, weights, boardSize, penalty] = getboard(tests.testsuite,whichBoards(n));

    time0 = cputime;
    board = solver(words, weights, boardSize, penalty);
    solverTime(n) = cputime-time0;

    [remainInDictionary, inPuzzle, bogus] = ...
        runSolution(board, boardSize, words, weights);

    solverScore(n) = scoresolution(remainInDictionary, bogus, penalty);

    % Print a report to the Command Window
    if flagVisualize
        visualizeSolution(board, boardSize, bogus, charFlag, fh);
        fprintf('\n  Puzzle: %d\n  Score: %d\n\n', n, solverScore(n));

        fprintf('Unused words remaining in dictionary\n------------------------------------\n');
        printWords(remainInDictionary.words, charFlag);

        fprintf('Valid words used in puzzle solution\n-----------------------------------\n');
        printWords(inPuzzle.words, charFlag);

        fprintf('Bogus words used in puzzle solution\n-----------------------------------\n');
        printWords(bogus.words, charFlag);

        if (n < numBoards) && mypause(fh)
            break; % stop
        end
    end
end

% Report results.
result = sum(solverScore);
resultTime = sum(solverTime);
message = sprintf('results: %.2f\ntime: %.2f', result, resultTime);

end

function printWords(words, charFlag)
%printWords: displays words in cell array to the Command Window

for n = 1:numel(words)
    if charFlag
        fprintf([char(words{n}),'\n']);
    else
        fprintf([num2str(words{n}),'\n']);
    end
end
fprintf('\n');
end

function [remainInDictionary, inPuzzle, bogus] = runSolution(board, boardSize, dictionary, weights)
%runSolution: parses the solution returned by the solver to determine the
% words it contains. Three structures are returned:
%
%  * remainInDictionary.words: words that remain in the dictionary
%  * remainInDictionary.points: points for each unused dictionary word
%  * remainInDictionary.hash: hash value used internally for searching
%
%  * inPuzzle.words: valid dictionary words used in the solution
%  * inPuzzle.points: points for each valid word
%
%  * bogus.words: words in the solution that aren't in the dictionary
%  * bogus.locations: linear indices of all letters in each word

% Initialize outputs
remainInDictionary = struct('words',{dictionary}, 'weights',{weights}, ...
    'hash', {cellfun(@applyHash, dictionary)});
inPuzzle = struct('words',{{}}, 'weights',{[]});
bogus = struct('words',{{}}, 'locations',{{}});

% Find words going across and down. dir = 1 for across, 2 for down.
for dir = 1:2
    for n = 1:boardSize
        % Get the current row / column to process
        if dir == 1
            curLine = board(n,:);
        else
            curLine = board(:,n).';
        end

        % Split current line of the board by locations of zeros
        zeroLoc = find(curLine == 0);
        idx1 = [1, zeroLoc+1];
        idx2 = [zeroLoc-1, boardSize];

        % Process the words on the current line
        for m = 1:numel(idx1)
            % Get the linear indices that make up the current word
            idxRange = idx1(m):idx2(m);
            idxConst = repmat(n, 1, numel(idxRange));
            if dir == 1
                linearWordIdx = sub2ind([boardSize, boardSize], idxConst, idxRange);
            else
                linearWordIdx = sub2ind([boardSize, boardSize], idxRange, idxConst);
            end

            % Current word to process
            curWord = board(linearWordIdx);

            if isscalar(curWord)
                % Get the indices of the adjacent squares - above,
                % below, left, and right.
                [r, c] = ind2sub([boardSize, boardSize], linearWordIdx);
                rows = [r-1, r+1, r,   r];
                cols = [c,   c,   c-1, c+1];
                outOfRange = (rows < 1) | (rows > boardSize) | ...
                    (cols < 1) | (cols > boardSize);
                rows(outOfRange) = [];
                cols(outOfRange) = [];
                adjacentIdx = sub2ind([boardSize, boardSize], rows, cols);

                % Treat a lone floating letter as invalid both across
                % and down. We assume that the dictionary never
                % contains any single letter words.
                if all(board(adjacentIdx) == 0)
                    bogus = addBogusWord(bogus, curWord, linearWordIdx);
                end
            elseif ~isempty(curWord) % Skip empty words - nothing to process
                % Determine if the word is in the dictionary
                dictionaryLoc = findWord(curWord, remainInDictionary.words, ...
                    remainInDictionary.hash);
                if isempty(dictionaryLoc)
                    % Bogus word - add to bogus word list
                    bogus = addBogusWord(bogus, curWord, linearWordIdx);
                else
                    % Legit word - move to inPuzzle list
                    [remainInDictionary, inPuzzle] = moveValidWord( ...
                        remainInDictionary, inPuzzle, dictionaryLoc);
                end
            end
        end
    end
end
end

function val = applyHash(word)
% applyHash: A rudimentary hashing function to speed up looking for words
% in the dictionary

n = min(5, ceil(numel(word)/2));
m1 = [1e10, 1e9, 1e8, 1e7, 1e6];
m2 = [1e4, 1e3, 1e2, 1e1, 1e0];

val = 1e12*sum(word) + sum(word(1:n).*m1(1:n)) + sum(word(end-n+1:end).*m2(1:n));
end

function idx = findWord(word, dictionary, dictionaryHash)
% findWord: Find word in dictionary using hash table. An empty matrix is
% return if the word does not exist in the dictionary

% Initialize output - assume word won't be found
idx = [];

% Find possible matches
possibleIdx = find(dictionaryHash == applyHash(word));

% Check all possible matches. If none is an actual match, we'll return empty
for n = 1:numel(possibleIdx)
    if isequal(dictionary{possibleIdx(n)},word)
        idx = possibleIdx(n);
        return;
    end
end

end

function bogus = addBogusWord(bogus, newWord, wordLoc)
%addBogusWord: add a word to the list of bogus words

bogus.words{end+1} = newWord;
bogus.locations{end+1} = wordLoc;
end

function [inDictionary, inPuzzle] = moveValidWord(inDictionary, inPuzzle, dictionaryLoc)
%moveValidWord: moves a word from the list of unused dictionary words to
% the list of words used in the puzzle solution

% Add word to list of words in the puzzle solution
inPuzzle.words{end+1} = inDictionary.words{dictionaryLoc};
inPuzzle.weights(end+1) = inDictionary.weights(dictionaryLoc);

% Remove word for list of unused dictionary words
inDictionary.words(dictionaryLoc) = [];
inDictionary.weights(dictionaryLoc) = [];
inDictionary.hash(dictionaryLoc) = [];
end

function score = scoresolution(remainInDictionary, bogus, penalty)
% SCORESOLUTION Calculates the score for the solution.

score = sum(remainInDictionary.weights) + penalty * numel(bogus.words);
end

function  [words, weights, n, penalty] = getboard(testsuite,k)
% GETBOARD Gets the k-th problem from the loaded testsuite

k = round(min(max(1,k),numel(testsuite)));
n = testsuite(k).n;
words = testsuite(k).words;
weights = testsuite(k).weights;
penalty = testsuite(k).penalty;
end

function visualizeSolution(board, boardSize, bogus, charFlag, fh)
% Display the resulting crossword puzzle solution

% Set the figure
figure(fh);
clf(fh);

% Set the axes
ax = axes('Parent',fh, 'Ydir','reverse', 'XLim',[0.5, boardSize+0.5], ...
    'YLim',[0.5, boardSize+0.5], 'Visible', 'off', ...
    'DataAspectRatio',[1, 1, 1], 'DataAspectRatioMode','manual');

% Draw the grid
for n = 0.5:boardSize+0.5
    line([n, n], [0.5, boardSize+0.5], 'Parent',ax, 'Color',[0, 0, 0]);
    line([0.5, boardSize+0.5], [n, n], 'Parent',ax, 'Color',[0, 0, 0]);
end

% Mark bogus word squares in red
for n = 1:numel(bogus.locations)
    for m = 1:numel(bogus.locations{n})
        [y, x] = ind2sub([boardSize, boardSize], bogus.locations{n}(m));
        patch([x-0.5, x+0.5, x+0.5, x-0.5],[y-0.5, y-0.5, y+0.5, y+0.5], [-1, -1, -1, -1], ...
            'Parent',ax, 'FaceColor',[1, 0, 0]);
    end
end

% Fill in each square - either a letter or black if unused
idx = 0;
for n = 1:boardSize
    for m = 1:boardSize
        if board(n,m) == 0
        patch([m-0.5, m+0.5, m+0.5, m-0.5],[n-0.5, n-0.5, n+0.5, n+0.5], [0, 0, 0, 0], ...
            'Parent',ax, 'FaceColor',[0, 0, 0]);
        else
            if charFlag
                letter = char(board(n,m));
            else
                letter = num2str(board(n,m));
            end

            idx = idx + 1;
            letterSize(idx) = 0.75/(boardSize*ceil(numel(letter)/2)); %#ok<AGROW>
            tx(idx) = text(m, n, letter, 'Parent',ax, 'FontUnits','normalized', ...
                'FontSize',letterSize(idx), 'VerticalAlignment','middle', ...
                'HorizontalAlignment','center', 'Clipping','on'); %#ok<AGROW>
        end
    end
end

% Set up zoom callback to resize text objects
h = zoom(fh);
set(h,'ActionPostCallback',@(o,e)zoomcb);

    function zoomcb
        zoom_ratio = boardSize / diff(get(ax,'XLim'));
        set(tx,{'FontSize'},num2cell(letterSize*zoom_ratio).');
    end
end

function flag = mypause(fh)
%mypause: provide navigation through puzzle board solutions

flag = false;
h1 = uicontrol('Parent',fh, 'Position',[20, 5, 100, 20], 'String','Next Problem', ...
    'Callback',@(o,e)uiresume(fh));
h2 = uicontrol('Parent',fh, 'Position',[140, 5, 100, 20], 'String','Stop', ...
    'Callback',@(o,e)stop_cb);
uiwait(fh);

% Remove uicontrols when we're done
if ishandle(h1)
    delete(h1);
end
if ishandle(h2)
    flag = ~isempty(get(h2,'UserData')); % stop?
    delete(h2);
end

    function stop_cb
        set(h2,'UserData',1);
        uiresume(fh);
    end
end


text:
function bytesize(in, fid)
% BYTESIZE writes the memory usage of the provide variable to the given file
% identifier. Output is written to screen if fid is 1, empty or not provided.

if nargin == 1 || isempty(fid)
    fid = 1;
end

s = whos('in');
fprintf(fid,[Bytes2str(s.bytes) '\n']);
end

function str = Bytes2str(NumBytes)
% BYTES2STR Private function to take integer bytes and convert it to
% scale-appropriate size.

scale = floor(log(NumBytes)/log(1024));
switch scale
    case 0
        str = [sprintf('%.0f',NumBytes) ' b'];
    case 1
        str = [sprintf('%.2f',NumBytes/(1024)) ' kb'];
    case 2
        str = [sprintf('%.2f',NumBytes/(1024^2)) ' Mb'];
    case 3
        str = [sprintf('%.2f',NumBytes/(1024^3)) ' Gb'];
    case 4
        str = [sprintf('%.2f',NumBytes/(1024^4)) ' Tb'];
    case -inf
        % Size occasionally returned as zero (eg some Java objects).
        str = 'Not Available';
    otherwise
       str = 'Over a petabyte!!!';
end
end


text:
function DoFSERecon(Simuh)

global VSig;
global VCtl;
global VObj;
global VCoi;
global VImg;

SX=reshape(VSig.Sx,VCtl.ResFreq * VCtl.FSE_ETL,VCtl.FSE_ShotNum,VCtl.SliceNum,VCoi.RxCoilNum,VObj.TypeNum); % matlab col priority
SY=reshape(VSig.Sy,VCtl.ResFreq * VCtl.FSE_ETL,VCtl.FSE_ShotNum,VCtl.SliceNum,VCoi.RxCoilNum,VObj.TypeNum);

SX=sum(SX, 5); % sum signal from all spin types
SY=sum(SY, 5);
%% FSE k-space recon
% sort k-space & row flipping
Sx=zeros(VCtl.ResFreq,VCtl.FSE_ETL*VCtl.FSE_ShotNum,VCtl.SliceNum,VCoi.RxCoilNum);
Sy=zeros(VCtl.ResFreq,VCtl.FSE_ETL*VCtl.FSE_ShotNum,VCtl.SliceNum,VCoi.RxCoilNum);
for e = 1: VCtl.FSE_ETL
    Sx(:,(e-1)*VCtl.FSE_ShotNum+1: e*VCtl.FSE_ShotNum, :,:) = reshape(SX((e-1)*VCtl.ResFreq+1: e*VCtl.ResFreq, :,:,:), [VCtl.ResFreq ,VCtl.FSE_ShotNum ,VCtl.SliceNum,VCoi.RxCoilNum]);
    Sy(:,(e-1)*VCtl.FSE_ShotNum+1: e*VCtl.FSE_ShotNum, :,:) = reshape(SY((e-1)*VCtl.ResFreq+1: e*VCtl.ResFreq, :,:,:), [VCtl.ResFreq ,VCtl.FSE_ShotNum ,VCtl.SliceNum,VCoi.RxCoilNum]);
end

% default -KxMax -> KxMax, -KyMax -> KyMax
% remove the most positive Kx point for making even number of Kx sample points (used for Matlab default fft)
Sx(end,:,:,:)=[];
Sy(end,:,:,:)=[];
ResFreq = VCtl.ResFreq - 1;

if isfield(VCtl,'ZF_Kx') % Zero Filling of k-space for increasing matrix size, no increae of resolution
    Sx2=zeros(str2double(VCtl.ZF_Kx),str2double(VCtl.ZF_Ky),str2double(VCtl.ZF_Kz(2:end))*VCtl.SliceNum,VCoi.RxCoilNum);
    Sy2=zeros(str2double(VCtl.ZF_Kx),str2double(VCtl.ZF_Ky),str2double(VCtl.ZF_Kz(2:end))*VCtl.SliceNum,VCoi.RxCoilNum);
    for i = 1:VCoi.RxCoilNum
        Sx2(str2double(VCtl.ZF_Kx)/2-ResFreq/2+1:str2double(VCtl.ZF_Kx)/2-ResFreq/2+1+ResFreq-1, ...
            str2double(VCtl.ZF_Ky)/2-(VCtl.FSE_ETL*VCtl.FSE_ShotNum)/2+1:str2double(VCtl.ZF_Ky)/2-(VCtl.FSE_ETL*VCtl.FSE_ShotNum)/2+1+(VCtl.FSE_ETL*VCtl.FSE_ShotNum)-1, ...
            str2double(VCtl.ZF_Kz(2:end))*floor(VCtl.SliceNum/2)-floor(VCtl.SliceNum/2)+1:str2double(VCtl.ZF_Kz(2:end))*floor(VCtl.SliceNum/2)-floor(VCtl.SliceNum/2)+1+VCtl.SliceNum-1,i)=Sx(:,:,:,i);
        Sy2(str2double(VCtl.ZF_Kx)/2-ResFreq/2+1:str2double(VCtl.ZF_Kx)/2-ResFreq/2+1+ResFreq-1, ...
            str2double(VCtl.ZF_Ky)/2-(VCtl.FSE_ETL*VCtl.FSE_ShotNum)/2+1:str2double(VCtl.ZF_Ky)/2-(VCtl.FSE_ETL*VCtl.FSE_ShotNum)/2+1+(VCtl.FSE_ETL*VCtl.FSE_ShotNum)-1, ...
            str2double(VCtl.ZF_Kz(2:end))*floor(VCtl.SliceNum/2)-floor(VCtl.SliceNum/2)+1:str2double(VCtl.ZF_Kz(2:end))*floor(VCtl.SliceNum/2)-floor(VCtl.SliceNum/2)+1+VCtl.SliceNum-1,i)=Sy(:,:,:,i);
    end
    Sx=Sx2;
    Sy=Sy2;
    clear Sx2 Sy2;
end
Sx=permute(Sx,[2 1 3 4]);
Sy=permute(Sy,[2 1 3 4]);

% match XYZ orientation
switch VCtl.ScanPlane
    case 'Axial'
        if  strcmp(VCtl.FreqDir,'A/P')
            Sx=permute(Sx,[2 1 3 4]);
            Sy=permute(Sy,[2 1 3 4]);
        else
            Sx=permute(Sx,[1 2 3 4]);
            Sy=permute(Sy,[1 2 3 4]);
        end
    case 'Sagittal'
        if  strcmp(VCtl.FreqDir,'S/I')
            Sx=permute(Sx,[2 1 3 4]);
            Sy=permute(Sy,[2 1 3 4]);
        else
            Sx=permute(Sx,[1 2 3 4]);
            Sy=permute(Sy,[1 2 3 4]);
        end
    case 'Coronal'
        if  strcmp(VCtl.FreqDir,'L/R')
            Sx=permute(Sx,[1 2 3 4]);
            Sy=permute(Sy,[1 2 3 4]);
        else
            Sx=permute(Sx,[2 1 3 4]);
            Sy=permute(Sy,[2 1 3 4]);
        end
end

S=Sx+1i*Sy;
for i = 1: VCoi.RxCoilNum
    VImg.Real(:,:,:,i)=real(fftshift(ifftn(fftshift(S(:,:,:,i)))));
    VImg.Imag(:,:,:,i)=imag(fftshift(ifftn(fftshift(S(:,:,:,i)))));
    VImg.Mag(:,:,:,i)=abs(VImg.Real(:,:,:,i)+1i*VImg.Imag(:,:,:,i));
    VImg.Phase(:,:,:,i)=angle(VImg.Real(:,:,:,i)+1i*VImg.Imag(:,:,:,i));
end
VImg.Sx(:,:,:,:)=Sx;
VImg.Sy(:,:,:,:)=Sy;

guidata(Simuh.SimuPanel_figure, Simuh);

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function [points3d, errors] = EsimatePosAmers(pointTracks, ...
    camPoses, cameraParams)

numTracks = numel(pointTracks);
points3d = zeros(numTracks, 3);
numCameras = size(camPoses, 2);
cameraMatrices = containers.Map('KeyType', 'uint32', 'ValueType', 'any');
for i = 1:numCameras
    id = camPoses(i).ViewId;
    R  = camPoses(i).Orientation;
    t  = camPoses(i).Location;
    size_t = size(t);
    if size_t(1) == 3
        t = t';
    end
    cameraMatrices(id) = cameraMatrix(cameraParams, R', -t*R');
end

for i = 1:numTracks
    track = pointTracks(i);
    points3d(i, :) = triangulateOnePoint(track, cameraMatrices);
end

if nargout > 1
    [~, errors] = reprojectionErrors(points3d, cameraMatrices, pointTracks);
end

%--------------------------------------------------------------------------
function point3d = triangulateOnePoint(track, cameraMatrices)

% do the triangulation
numViews = numel(track.ViewIds);
A = zeros(numViews * 2, 4);
for i = 1:numViews
    id = track.ViewIds(i);
    P = cameraMatrices(id)';
    A(2*i - 1, :) = track.Points(i, 1) * P(3,:) - P(1,:);
    A(2*i    , :) = track.Points(i, 2) * P(3,:) - P(2,:);
end

[~,~,V] = svd(A);
X = V(:, end);
X = X/X(end);
point3d = X(1:3)';

%--------------------------------------------------------------------------
function [errors, meanErrorsPerTrack] = reprojectionErrors(points3d, ...
    cameraMatrices, tracks)
numPoints = size(points3d, 1);
points3dh = [points3d, ones(numPoints, 1)];
meanErrorsPerTrack = zeros(numPoints, 1);
errors = [];
for i = 1:numPoints
    p3d = points3dh(i, :);
    reprojPoints2d = reprojectPoint(p3d, tracks(i).ViewIds, cameraMatrices);
    e = sqrt(sum((tracks(i).Points - reprojPoints2d).^2, 2));
    meanErrorsPerTrack(i) = mean(e);
    errors = [errors; e];
end

%--------------------------------------------------------------------------
function points2d = reprojectPoint(p3dh, viewIds, cameraMatrices)
numPoints = numel(viewIds);
points2d = zeros(numPoints, 2);
for i = 1:numPoints
    p2dh = p3dh * cameraMatrices(viewIds(i));
    points2d(i, :) = p2dh(1:2) ./ p2dh(3);
end

author: mbrossar


text:
classdef Voigt2TensorPSrepresentation < handle

    properties (Access = protected)
        indexTransformer
    end

    methods (Access = protected)

        function createIndexTransformer(obj)
            obj.indexTransformer = TensorVoigtIndexTransformerPS();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function result = isInGeodesicROI(mesh_total,wholepath,face)

wholepath=wholepath(2:end);
verticesinfaces=mesh_total.faces(face,:);
b = intersect(verticesinfaces,wholepath);

resultat = size(b,2) < 1;
result = resultat(1);


text:
function irepistruct = irepifit(irepistruct)
% irepistruct = irepifit(irepistruct)
% Performs LMS fit.
%

s = irepistruct; % copy into s for easy handling

% s = irepitiming(s); Need to do this if fitting seq params

% Need to be able to cache so not recomputing every slice
% Need this way too for incremental
s = irepisynth(s);

s.X0 = s.yRead;

s.iexclude1 = [];
if(s.nexclude > 0)
  % Just skip the first n
  % s.iexclude1 = [1:s.nexclude]';

  % This skips time points according to a rule. If a given time
  % point was the last acquired slice, then the next time point
  % for that slice is skipped.
  %indSlice = find((s.EventSliceNo == s.sliceno | s.EventSliceNo < 0) & s.IsReadOut);
  %nthSliceAcq = s.AcqSliceNo(indSlice);
  %indLast = find(nthSliceAcq == s.nslices);
  %if(~isempty(indLast))
  %  indEx = indLast + 1;
  %  ok = find(indEx <= s.ntp);
  %  indEx = indEx(ok);
  %else
  %  indEx = [];
  %end
  %s.iexclude1 = indEx;

  % This skips time points of the first acquired slice after inversion
  indSlice = find(s.EventSliceNo == s.sliceno | s.EventSliceNo < 0);
  indROS = find(s.IsReadOut(indSlice));
  nthSliceAcq = s.AcqSliceNo(indSlice(indROS));
  s.iexclude1 = find(nthSliceAcq <= s.nexclude);
  %s.iexclude1 = find(nthSliceAcq > s.nexclude);
end
s.iexclude2 = [];
if(s.nminexclude > 0)
  [ysort isort] = sort(s.y);
  s.iexclude2 = sort(isort(1:s.nminexclude));
end

if(~isempty(s.iexclude1) & isempty(s.iexclude2))
  s.iexclude = s.iexclude1;
elseif(isempty(s.iexclude1) & ~isempty(s.iexclude2))
  s.iexclude = s.iexclude2;
else
  s.iexclude = unique([s.iexclude1; s.iexclude2]);
end
s.indkeep = setdiff([1:s.ntp],s.iexclude);

% remove the mean
%s.X0 = s.X0-repmat(mean(s.X0),[size(s.X0,1) 1]);
%s.y = s.y-mean(s.y);

s.yFit = s.y(s.indkeep,:);
s.X = s.X0(s.indkeep,:);

s.tFit = s.tRead(s.indkeep);

nX = size(s.X,2);
Xvar = sum(s.X.*s.X); % need to cache, unless exclude
s.M0 = sum(s.X.*repmat(s.yFit,[1 nX]))./Xvar;
s.yhat = s.X.*repmat(s.M0,[length(s.yFit) 1]);
s.res = repmat(s.yFit,[1 nX])-s.yhat;
s.yhat0 = s.X0.*repmat(s.M0,[length(s.y) 1]);
if(0)
  s.M0 = inv(s.X'*s.X)*(s.X'*s.y);
  s.yhat = s.M0*s.X;
  s.res = s.y-s.yhat;
  s.yhat0 = s.M0*s.X0;
end
s.rstd = std(s.res);
err = s.rstd;
s.err = err;

irepistruct = s;

return;


text:
function cin_test ( )

%*****************************************************************************80
%
%% CIN_TEST tests R4_CIN and R8_CIN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CIN_TEST:\n' );
  fprintf ( 1, '  Test CIN_VALUES, R4_CIN, R8_CIN.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             X          CIN(X)\n' );
  fprintf ( 1, '                     R4_CIN(X)        Diff\n' );
  fprintf ( 1, '                     R8_CIN(X)        Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx1 ] = cin_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_cin ( single ( x ) );
    fx3 = r8_cin ( x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.6g\n', x, fx1 );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer


text:
function [rhsE, rhsH] = MaxwellRHS1D(E,H,eps,mu)

% function [rhsE, rhsH] = MaxwellRHS1D(E,H,eps,mu)
% Purpose  : Evaluate RHS flux in 1D Maxwell

Globals1D;

% Compute impedance
Zimp = sqrt(mu./eps);

% Define field differences at faces
dE = zeros(Nfp*Nfaces,K); dE(:) = E(vmapM)-E(vmapP);
dH = zeros(Nfp*Nfaces,K); dH(:) = H(vmapM)-H(vmapP);
Zimpm = zeros(Nfp*Nfaces,K); Zimpm(:) = Zimp(vmapM);
Zimpp = zeros(Nfp*Nfaces,K); Zimpp(:) = Zimp(vmapP);
Yimpm = zeros(Nfp*Nfaces,K); Yimpm(:) = 1./Zimpm(:);
Yimpp = zeros(Nfp*Nfaces,K); Yimpp(:) = 1./Zimpp(:);

% Homogeneous boundary conditions, Ez=0
Ebc = -E(vmapB); dE (mapB) = E(vmapB) - Ebc;
Hbc =  H(vmapB); dH (mapB) = H(vmapB) - Hbc;

% evaluate upwind fluxes
fluxE = 1./(Zimpm + Zimpp).*(nx.*Zimpp.*dH - dE);
fluxH = 1./(Yimpm + Yimpp).*(nx.*Yimpp.*dE - dH);

% compute right hand sides of the PDE's
rhsE = (-rx.*(Dr*H) + LIFT*(Fscale.*fluxE))./eps;
rhsH = (-rx.*(Dr*E) + LIFT*(Fscale.*fluxH))./mu;
return