text:
function van_der_corput_test04 ( )

%*****************************************************************************80
%
%% VAN_DER_CORPUT_TEST04 tests VAN_DER_CORPUT_SEQUENCE, VAN_DER_CORPUT_SEED_GET, VAN_DER_CORPUT_SEED_SET.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'VAN_DER_CORPUT_TEST04\n' );
  fprintf ( 1, '  VAN_DER_CORPUT_BASE_SET sets the base;\n' );
  fprintf ( 1, '  VAN_DER_CORPUT_BASE_GET gets the base;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  In this test, we show how setting the base changes\n' );
  fprintf ( 1, '  the results of the sequence.\n' );
  fprintf ( 1, '\n' );

  seed = 0;
  van_der_corput_seed_set ( seed );

  base = 2;
  van_der_corput_base_set ( base );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Current SEED = %d\n', seed );
  fprintf ( 1, '  Current BASE = %d\n', base );

  n = 10;
  r = van_der_corput_sequence ( n );

  for i = 1: n
    fprintf ( 1, '%d %f\n', seed+i-1, r(i) );
  end

  seed = 0;
  van_der_corput_seed_set ( seed );

  base = 3;
  van_der_corput_base_set ( base );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Repeat the calculation with a new base.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Current SEED = %d\n', seed );
  fprintf ( 1, '  Current BASE = %d\n', base );

  seed = van_der_corput_seed_get ( );

  n = 10;
  r = van_der_corput_sequence ( n );

  for i = 1: n
    fprintf ( 1, '%d %f\n', seed+i-1, r(i) );
  end

  seed = 0;
  van_der_corput_seed_set ( seed );

  base = 4;
  van_der_corput_base_set ( base );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Repeat the calculation with a new base.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Current SEED = %d\n', seed );
  fprintf ( 1, '  Current BASE = %d\n', base );

  seed = van_der_corput_seed_get ( );

  n = 10;
  r = van_der_corput_sequence ( n );

  for i = 1: n
    fprintf ( 1, '%d %f\n', seed+i-1, r(i) );
  end

  return
end


text:
function cCell = cqtSparse2Cell(cSparse,M,cDC,cNyq)

bins = size(M,1)/2 - 1;
cCell = cell(1,bins+2);
cCell{bins+2} = cNyq;

M = M(1:bins+1);
step = 1;
cSparse = full(cSparse);
distinctHops = log2(M(bins+1)/M(2))+1;
curNumCoef = M(bins+1);

for ii=1:distinctHops
   idx = (M == curNumCoef);
   temp = cSparse(idx,1:step:end).';
   temp = num2cell(temp,1);
   cCell(idx) = temp;
   step = step*2;
   curNumCoef = curNumCoef / 2;
end

cCell{1} = cDC;

author: nesl


text:
function DrawClusters(h1,M,gIX,dataFR,numK,stim,fictive,clrmap,rankscore,iswrite)
pos = get(gca,'Position');
barratio = 0.03;

%% Prepare cluster data
% down-sample
displaymax = 1000;
numcell = size(M,1);
if numcell > displaymax,
    skip = round(numcell/displaymax);
    M = M(1:skip:end,:);
    gIX = gIX(1:skip:end,:);
end

numK = double(max(double(numK),double(max(gIX))));

% sort traces by index
im = M;
[nLines,nFrames] = size(im);
[~,I] = sort(gIX);
im = im(I,:);

%% convert imagesc effect into rgb matrix 'RGB'
if dataFR, % is drawing cell-response fluorescence
    cmap = gray(64);
    im = AutoScaleImage0to1(im); % scaling to min/max 0/1
    minlim = 0;
    maxlim = 1;
else % create blue-white-red map, for regression results
    cmap = zeros(64,3);
    cmap(:,1) = [linspace(0,1,32), linspace(1,1,32)];
    cmap(:,2) = [linspace(0,1,32), linspace(1,0,32)];
    cmap(:,3) = [linspace(1,1,32), linspace(1,0,32)];
    minlim = -1; %min(min(im));
    maxlim = 1; %max(max(im));
end
RGB = ImageToRGB(im,cmap,minlim,maxlim); % map image matrix to range of colormap

%% add vertical color code
if exist('clrmap','var'),
    if strcmp(clrmap,'jet'),
        temp = flipud(jet(numK));
    else % 'hsv'
        temp = hsv(round(numK*1.1));
    end
else % 'hsv'
    temp = hsv(round(numK*1.1));
end
cmap2 = vertcat(temp(1:numK,:),[0,0,0]); % extend colormap to include black
bwidth = max(round(nFrames/30),1);

idx = gIX(I);
ix_div = [find(diff(idx));length(idx)];

bars = ones(nLines,bwidth);
if numK>1,
    bars(1:ix_div(1),:) = idx(ix_div(1));
    for i = 2:length(ix_div),
        % paint color
        bars(ix_div(i-1)+1:ix_div(i),:) = idx(ix_div(i));
    end
end
im_bars = reshape(cmap2(bars,:),[size(bars),3]);
% add a white division margin
div = ones(nLines,round(bwidth/2),3);
% put 3 parts together
im = horzcat(RGB,div,im_bars);

%% plot figure

[s1,s2,s3] = size(im);

hold off;
set(h1,'Position',[pos(1),pos(2)+pos(4)*barratio*3,pos(3),pos(4)*(1-4*barratio)]);
if s1<30,
    temp = im;
    im = ones(30,s2,s3);
    im(1:s1,:,:) = temp;
end
image(im);
set(gca, 'box', 'off')

hold on;

% plot cluster division lines
plot([0.5,s2+0.5],[0.5,0.5],'k','Linewidth',0.5);
if numK>1,
    for i = 1:length(ix_div),% = numK-1,
        y = ix_div(i)+0.5;
        plot([0.5,s2+0.5],[y,y],'k','Linewidth',0.5);
    end
end

ylabel(['Number of cells: ' num2str(numcell)]);
set(gca,'YTick',[],'XTick',[]);
set(gcf,'color',[1 1 1]);

colormap gray;

% write in number label
x = s2*1.003;
y_last = 0;
for i = 1:length(ix_div),
    % avoid label crowding
    margin = 0.015*s1;
    y0 = ix_div(i)+0.5;
    y = max(y_last+margin,y0);
    if i<length(ix_div),
        ynext = ix_div(i+1);
    else
        ynext = y0+margin*2;
    end
    % draw if not squeezed too much away
    if y < y0+margin*2 && y < 0.5*(y0+ynext), %nLines+stimheight-margin,
        if exist('iswrite','var') && iswrite,
            try
                text(x,y,[num2str(idx(ix_div(i))) ': ' num2str(rankscore(i))],'HorizontalAlignment','Left','VerticalAlignment','Bottom',...
                    'FontUnits','normalized','FontSize',0.015,'Tag','label');
            catch
                disp('rankscore invalid');
                text(x,y,num2str(idx(ix_div(i))),'HorizontalAlignment','Left','VerticalAlignment','Bottom',...
                    'FontUnits','normalized','FontSize',0.015,'Tag','label');
            end

        else
            text(x,y,num2str(idx(ix_div(i))),'HorizontalAlignment','Left','VerticalAlignment','Bottom',...
                'FontUnits','normalized','FontSize',0.015,'Tag','label');
        end
        y_last = y;
    end
end

%% Draw stim bars
halfbarheight = 100;
stimbar_ = GetStimBar(halfbarheight,stim);

% pad with margin on right to fit combo-image size
barlength = size(im,2);
stimbar = ones(halfbarheight*2,barlength,3);
stimbar(:,1:nFrames,:) = stimbar_(:,1:nFrames,:); % crop to size

% top axis
axes('Position',[pos(1),pos(2)+pos(4)*(1-barratio),pos(3),pos(4)*barratio]);
image(stimbar);axis off;hold on;
plot([1,length(stim)],[1,1],'k','Linewidth',0.5); % plot top border
plot([0.5,0.5],[0,size(stimbar,1)],'k','Linewidth',0.5); % plot left/right border
plot([length(stim),length(stim)],[0,size(stimbar,1)],'k','Linewidth',0.5);

% bottom axis
axes('Position',[pos(1),pos(2)+pos(4)*barratio*2,pos(3),pos(4)*barratio]);
image(stimbar);axis off;hold on;
plot([1,length(stim)],[size(stimbar,1),size(stimbar,1)],'k','Linewidth',0.5); % plot bottom border
plot([0.5,0.5],[0,size(stimbar,1)],'k','Linewidth',0.5); % plot left/right border
plot([length(stim),length(stim)],[0,size(stimbar,1)],'k','Linewidth',0.5);

%% Draw fictive
axes('Position',[pos(1),pos(2),pos(3),pos(4)*barratio*2]);

temp = ones(size(fictive,1),barlength);
temp(:,1:length(fictive)) = fictive;
imagesc(temp);colormap hot
set(gca,'YTick',[],'XTick',[]);
set(gcf,'color',[1 1 1]);
set(gca, 'box', 'off')
hold on;axis ij;

% plot division lines
for i = 0:3,
    y = i+0.5;
    plot([0.5,length(fictive)+0.5],[y,y],'w','Linewidth',0.5);
end
% labels
names = {'Left','Right','Forward','Raw L','Raw R'};
x = -s2*0.05;
for i = 1:5,
    y = i;
    text(x,y,names{i},'Fontsize',7);
end

%% horizontal axis: approximate time length, assuming 2 fps (real is 1.97 for now)
set(gca,'XColor',[1 1 1]);
if s2/2/60<5, % <5 min
    l = xlabel(['Length of display: ~ ' num2str(s2/2) 'sec'],'Color',[0 0 0]);
else
    l = xlabel(['Length of display: ~ ' num2str(s2/2/60) 'min'],'Color',[0 0 0]);
end

%% Draw fish icon:
if exist('fishpic.jpg','file')==2,
    fishpic = imread('fishpic.jpg');
    axes('Position',[0.015, 0.095, 0.03, 0.02]);
    image(fishpic); axis off; axis equal
end
end

%% Internal Functions

function im = AutoScaleImage0to1(im)
temp = mat2gray(im);
low_high = stretchlim(temp,[0.005,0.995]);
temp(temp<low_high(1)) = low_high(1);
temp(temp>low_high(2)) = low_high(2);
im = mat2gray(temp);
end

function RGB = ImageToRGB(im,cmap,minlim,maxlim)
L = size(cmap,1);
ix = round(interp1(linspace(minlim,maxlim,L),1:L,im,'linear','extrap'));
RGB = reshape(cmap(ix,:),[size(ix) 3]); % Make RGB image from scaled.
end

author: xiuyechen

repo: FishExplorer


text:
function classifier = train_boosted_dt_mc(features, cat_features, labels, ...
    num_iterations, num_nodes, stopval, init_weights, varargin)
%
%classifier = train_boosted_dt_mc(features, cat_features, labels, ...
%    num_iterations, num_nodes, stopval, init_weights, varargin)
%
% Train a classifier based on boosted decision trees.  Boosting done by the
% logistic regression version of Adaboost (Adaboost.L - Collins, Schapire,
% Singer 2002).  At each
% iteration, a set of decision trees is created for each class, with
% confidences equal to 1/2*ln(P+/P-) for that class, according to the
% weighted distribution.  Final classification is based on the largest
% confidence label (possibly incorporating a prior as h0(c) =
% 1/2*ln(Pc/(1-Pc)).  Weights are assigned as
% w(i,j) = 1 / (1+exp(sum{t in iterations}[yij*ht(xi, j)])).

if length(varargin) == 1  % class names supplied
    gn = varargin{1};
    gid = zeros(size(labels));
    for c = 1:length(gn)
        ind = find(strcmp(labels, gn{c}));
        gid(ind) = c;
        if ~isempty(init_weights)
            disp([gn{c} ': ' num2str(sum(init_weights(ind)))]);
        else
            disp([gn{c} ': ' num2str(length(ind))]);
        end
    end
    ind = find(gid==0);
    gid(ind) = [];
    labels(ind) = [];
    features(ind, :) = [];
else
    [gid, gn] = grp2idx(labels);
end

if ~exist('stopval', 'var') || isempty(stopval)
    stopval = 0;
end
if ~exist('init_weights', 'var')
    init_weights = [];
end

classifier.names = gn;

num_classes = length(gn);
num_data = length(gid);

if isempty(init_weights)
    init_weights = ones(num_data, 1)/num_data;
else
    init_weights = init_weights / sum(init_weights);
end

% if no examples from a class are present, create one dummy example for
% that class with very small weight
for c = 1:numel(gn)
    if ~any(gid==c)
        disp(['warning: no examples from class ' gn(c)])
        gid(end+1) = c;
        features(end+1, :) = zeros(size(features(end, 1)));
        num_data = num_data + 1;
        init_weights(end+1) = min(init_weights)/2;
    end
end

all_conf = zeros(num_data, num_classes);
for c = 1:num_classes

    disp(['class: ' num2str(gn{c})]);
    y = (gid == c)*2-1;
    cl = [-1 1];
    nc = 2;
    w = zeros(num_data, 1);
    cw = zeros(num_classes, 1);
    for i = 1:2
        indices = find(y==cl(i));
        %count = sum(init_weights(indices));
        %w(indices) = init_weights(indices) / count / 2;
        w(indices) = init_weights(indices);

        if cl(i)==1
            %classifier.h0(c) = log(count / (1-count));
            classifier.h0(c) = 0;
        end

    end

    data_confidences = zeros(num_data, 1);
    aveconf = [];

    for t = 1:num_iterations
        % learn decision tree based on weighted distribution
        dt = treefitw(features, y, w, 1/num_data/2, 'catidx', cat_features, 'method', 'classification', 'maxnodes', num_nodes*4);
        [tmp, level] = min(abs(dt.ntermnodes-num_nodes));
        dt = treeprune(dt, 'level', level-1);

        % assign partition confidences
        pi = (strcmp(dt.classname{1},'1')) + (2*strcmp(dt.classname{2},'1'));
        ni = (strcmp(dt.classname{1},'-1')) + (2*strcmp(dt.classname{2},'-1'));

        classprob = dt.classprob;
        confidences = 1/2*(log(classprob(:, pi)) - log(classprob(:, ni)));

        % assign weights
        [class_indices, nodes, classes] = treeval(dt, features);
        data_confidences = data_confidences + confidences(nodes);

        w = 1 ./ (1+exp(y.*data_confidences));
        w = w / sum(w);

        %disp(['c: ' num2str(mean(1 ./ (1+exp(-y.*data_confidences)))) '  e: ' num2str(mean(y.*data_confidences < 0)) '   w: ' num2str(max(w))]);

        classifier.wcs(t, c).dt = dt;
        classifier.wcs(t, c).confidences = confidences;

        aveconf(t) = mean(1 ./ (1+exp(-y.*data_confidences)));
        if t>10 && (aveconf(t)-aveconf(t-10) < stopval)
            disp(num2str(aveconf))
            disp(['Stopping after ' num2str(t) ' trees'])
            break;
        end

    end

    finalconf = 1 ./ (1+exp(-y.*data_confidences));
    finalerr = (y.*data_confidences < 0);
    disp(['confidence:: mean: ' num2str(mean(finalconf)) ...
        '  pos: ' num2str(mean(finalconf(y==1))) ...
        '  neg: ' num2str(mean(finalconf(y~=1)))]);
    disp(['training error:: mean: ' num2str(mean(finalerr)) ...
        '  pos: ' num2str(mean(finalerr(y==1))) ...
        '  neg: ' num2str(mean(finalerr(y~=1)))]);
    all_conf(:, c) = data_confidences+classifier.h0(c);

end

% compute and display training error
[tmp, assigned_label] = max(all_conf, [], 2);
conf_matrix = zeros(num_classes, num_classes);
for c = 1:num_classes
    indices = find(gid==c);
    for c2 = 1:num_classes
        conf_matrix(c, c2) = mean(assigned_label(indices)==c2);
    end
    disp([gn{c} ' error: ' num2str(mean(assigned_label(indices)~=c))]);
end
disp('Confusion Matrix: ');
disp(num2str(conf_matrix));
disp(['total error: ' num2str(mean(assigned_label~=gid))]);

author: zouchuhang

repo: LayoutNet


text:
function out=first_order_reshape(S,ref)

S=S{2};

for j=1:size(S.signal,2)
out(j)= S.signal{j} / S.signal{ref};
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
%% This function will calculate the prediction error between the discovered system and actual test data
% There are two methods to calculate the prediction error.
%
% Method one: This method suits for the one dimensional system only. The idea is
% using ODE45 to drive each point k-steps forward in time and calculate the
% prediction error between the system estimation and actual states. Here,
% the system's estimate is the value of derivative at k-steps.
%
% Method two: This method suits for one dimensional system and
% multi-dimensional system. We use immediate prediction of the system's
% derivative and calculate the prediction error of the system model's
% derivative output and actual derivative.
%
% Last Updated: 2019/06/04
% Coded By: K
%%
function [Score]=Get_Score(dData_test,Data_test,u,Control,tspan,Shuffle,name,Prediction_Steps,dt,size1,size2,method)
% Get the ODE simulation result
Noise=0;

% Start calculate
if method ==1
    % If the method one is selected, then we use the multi step prediction

    % Generate the function handel
    if u==0
        ODE_func=str2func(strcat('@(t,z)',name,'(t,z)'));
    end

    % Pass the test data as dummy variable
    Dummy=Data_test;

    if Prediction_Steps==0
        % This will result immediate prediction
        dData_Es=ODE_func(0,Dummy);
        Error=norm(dData_test-dData_Es);
        Dem=norm(dData_test);
        Score=Error/Dem;
    else
        dData_Es=ODE_func(0,Dummy);
        % Get the function estimation using RK45
        for pinpin=1:Prediction_Steps
            RK_k1=dt*dData_Es;
            RK_k2=dt*ODE_func(0,Dummy+0.5*RK_k1);
            RK_k3=dt*ODE_func(0,Dummy+0.5*RK_k2);
            RK_k4=dt*ODE_func(0,Dummy+RK_k3);
            Dummy=Dummy+(1/6)*(RK_k1+2*RK_k2+2*RK_k3+RK_k4);
            dData_Es=ODE_func(0,Dummy);
        end

        % Reshape
        Dum1=reshape(dData_Es,size1,size2);
        Dum2=reshape(dData_test,size1,size2);

        % Arrange all the data to the same time
        Dum3=Dum1(1:end-Prediction_Steps,:);
        Dum4=Dum2(1+Prediction_Steps:end,:);

        % Calculate the prediction error
        Error=norm(reshape((Dum4-Dum3),[],1));
        Dem=norm(reshape(Dum4,[],1));

        % Calculate the norm
        Score=Error/Dem;
    end
else
    % Generate the function handel
    if u==0
        ODE_func=str2func(strcat('@(t,z)',name,'(t,z)'));
    end

    % Get the function estimation
    dData_Es=ODE_func(0,Data_test);
    Error=norm(dData_test-dData_Es);
    Dem=norm(dData_test);
    Score=Error/Dem;
end

author: dynamicslab

repo: SINDy-PI


text:
function res = demo_DCFNet()
addpath(fullfile('..','DCFNet'));

init_rect = [129,80,64,78];
img_file = dir('./David/img/*.jpg');
img_file = fullfile('./David/img/', {img_file.name});
subS.init_rect = init_rect;
subS.s_frames = img_file;

param = [];
param.gpu = true;
gpuDevice(1);
param.visual = true;

res = run_DCFNet(subS,0,0,param);
disp(['fps: ',res.fps]);

end

author: foolwood

repo: DCFNet

save_path: github-repos/MATLAB/foolwood-DCFNet


text:
function [HDR]=scpopen(arg1,CHAN,arg4,arg5,arg6)
% SCPOPEN reads and writes SCP-ECG files
%
% SCPOPEN is an auxillary function to SOPEN for
% opening of SCP-ECG files for reading ECG waveform data
%
% Use SOPEN instead of SCPOPEN
%
% See also: fopen, SOPEN,

%	$Id: scpopen.m 2205 2009-10-27 12:18:15Z schloegl $
%	(C) 2004,2006,2007,2008 by Alois Schloegl <a.schloegl@ieee.org>
%    	This is part of the BioSig-toolbox http://biosig.sf.net/
%
%    BioSig 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.

if nargin<2, CHAN=0; end;

if isstruct(arg1)
        HDR=arg1;
        FILENAME=HDR.FileName;
elseif ischar(arg1);
        HDR.FileName=arg1;
        fprintf(2,'Warning SCPOPEN: the use of SCPOPEN is discouraged; please use SOPEN instead.\n');
end;

VER = version;

fid = fopen(HDR.FileName,HDR.FILE.PERMISSION,'ieee-le');
HDR.FILE.FID = fid;
if ~isempty(findstr(HDR.FILE.PERMISSION,'r')),		%%%%% READ
	tmpbytes = fread(fid,inf,'uchar');
        tmpcrc   = crc16eval(tmpbytes(3:end));
	fseek(fid, 0, 'bof');
        HDR.FILE.CRC = fread(fid,1,'uint16');
	if (HDR.FILE.CRC ~= tmpcrc);
		fprintf(HDR.FILE.stderr,'Warning: CRC check failed (%x vs %x)\n',tmpcrc,HDR.FILE.CRC);
        end;

        HDR.FILE.Length = fread(fid,1,'uint32');
        if HDR.FILE.Length~=HDR.FILE.size,
                fprintf(HDR.FILE.stderr,'Warning SCPOPEN: header information contains incorrect file size %i %i \n',HDR.FILE.Length,HDR.FILE.size);
        end;
	HDR.data = [];

        DHT = [0,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,7,-7,8,-8,9,-9;0,1,5,3,11,7,23,15,47,31,95,63,191,127,383,255,767,511,1023]';
        prefix  = [1,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,10,10];
	PrefixLength = prefix;

        %PREFIX = [0,4,5,12,13,28,29,60,61,124,125,252,253,508,509,1020,1021,1022,1023];
        PREFIX  = [0,4,5,12,13,28,29,60,61,124,125,252,253,508,509,1020,1021,1022,1023]'.*2.^[32-prefix]';
        codelength = [1,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,18,26];
        mask    = [1,7,7,15,15,31,31,63,63,127,127,255,255,511,511,1023,1023,1023,1023]'.*2.^[32-prefix]';
        %MASK    = dec2bin(mask);
        %mask   = [1,7,7,15,15,31,31,63,63,127,127,255,255,511,511,1023,1023,1023,1023]';

        mask2    = [1,7,7,15,15,31,31,63,63,127,127,255,255,511,511,1023,1023,1023,1023]';
        PREFIX2  = DHT(:,2);

	HT19999 = [prefix',codelength',ones(length(prefix),1),DHT];
	HT = [prefix',codelength',ones(length(prefix),1),DHT];

        dd = [0:255]';
        ACC = zeros(size(dd));
        c = 0;
        for k2 = 1:8,
                ACC = ACC + (dd>127).*(2^c);
                dd  = mod(dd*2, 256);
                c   = c + 1;
        end;

        section.CRC     = fread(fid,1,'uint16');
        section.ID      = fread(fid,1,'uint16');
        section.Length  = fread(fid,1,'uint32');
        section.Version = fread(fid,[1,2],'uint8');
        section.tmp     = fread(fid,[1,6],'uint8');

        NSections = min(11,(section.Length-16)/10);
        for k = 1:NSections,
                HDR.Block(k).id = k;
                HDR.Block(k).length = 0;
                HDR.Block(k).startpos = -1;
        end;
	for K = 1:NSections,
                k = fread(fid,1,'uint16');
		len = fread(fid,1,'uint32');
		pos = fread(fid,1,'uint32');
		if ((k > 0) && (k < NSections))
	                HDR.Block(k).id = k;
        	        HDR.Block(k).length = len;
                	HDR.Block(k).startpos = pos-1;

%% [HDR.Block(k).id ,length(tmpbytes), HDR.Block(k).length, HDR.Block(k).length+HDR.Block(k).startpos]
		%% FIXME: instead of min(...,FileSize) a warning or error message should be reported
	                tmpcrc = crc16eval(tmpbytes(HDR.Block(k).startpos+3:min(HDR.Block(k).startpos+HDR.Block(k).length,HDR.FILE.size)));

	                if (HDR.Block(k).length>0)
        	        if (tmpcrc~=(tmpbytes(HDR.Block(k).startpos+(1:2))'*[1;256]))
                		fprintf(HDR.FILE.stderr,'Warning SCPOPEN: faulty CRC %04x in section %i\n',tmpcrc,k-1);
	                end;
        	        end;
		end;
        end;

%%[[HDR.Block.id];[HDR.Block.length];[HDR.Block.startpos]]'

	% default values - in case Section 6 is missing
	HDR.NS = 0; HDR.SPR = 0; HDR.NRec = 0; HDR.Calib = zeros(1,0);
        secList = find([HDR.Block.length]);
        for K = secList(1:end),
                if fseek(fid,HDR.Block(K).startpos,'bof');
                        fprintf(HDR.FILE.stderr,'Warning SCPOPEN: section %i not available, although it is listed in Section 0\n',secList(K+1));
                end;
                section.CRC     = fread(fid,1,'uint16');
                section.ID      = fread(fid,1,'uint16');
                section.Length  = fread(fid,1,'uint32');
                section.Version = fread(fid,[1,2],'uint8');
                section.tmp     = fread(fid,[1,6],'uint8');

                HDR.SCP.Section{find(K==secList)} = section;
                if (section.Length==0),
                elseif section.ID==0,
                        NSections = (section.Length-16)/10;
                        for k = 1:NSections,
                                HDR.Block(k).id = fread(fid,1,'uint16');
                                HDR.Block(k).length = fread(fid,1,'uint32');
                                HDR.Block(k).startpos = fread(fid,1,'uint32')-1;
                        end;

                elseif section.ID==1,
                        tag = 0;
                        k1  = 0;
                        Sect1Len = section.Length-16;
                        ListOfRequiredTags = [2,14,25,26];
                        ListOfRecommendedTags = [0,1,5,8,15,34];
                        while (tag~=255) & (Sect1Len>2),
                                tag = fread(fid,1,'uint8');
                                len = fread(fid,1,'uint16');
                                Sect1Len = Sect1Len - 3 - len;
%% [tag,len,Sect1Len],          %% DEBUGGING information
                                field = fread(fid,[1,len],'uchar');
                                if tag == 0,
                                        ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
                                        HDR.Patient.Name = char(field);  %% LastName
                                elseif tag == 1,
                                        ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
                                        HDR.Patient.FirstName = char(field);
                                elseif tag == 2,
                                        ListOfRequiredTags(find(ListOfRequiredTags==2))=[];
                                        HDR.Patient.Id = char(field);
                                elseif tag == 3,
                                        HDR.Patient.LastName2 = char(field);
                                elseif tag == 4,
                                        HDR.Patient.Age = field(1:2)*[1;256];
                                        tmp = field(3);
                                        if     tmp==1, HDR.Patient.Age = HDR.Patient.Age; % unit='Y';
                                        elseif tmp==2, HDR.Patient.Age = HDR.Patient.Age/12; % unit='M';
                                        elseif tmp==3, HDR.Patient.Age = HDR.Patient.Age/52; % unit='W';
                                        elseif tmp==4, HDR.Patient.Age = HDR.Patient.Age/365.25; % unit='d';
                                        elseif tmp==5, HDR.Patient.Age = HDR.Patient.Age/(365.25*24); %unit='h';
                                        else warning('units of age not specified');
                                        end;
                                elseif (tag == 5)
                                        ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
                                	if any(field(1:4)~=0)
                                        	HDR.Patient.Birthday = [field(1:2)*[1;256],field(3:4),12,0,0];
                                        end;
                                elseif (tag == 6)
                                	if any(field(1:3)),
                                        HDR.Patient.Height = field(1:2)*[1;256];
                                        tmp = field(3);
                                        if tmp==1, % unit='cm';
                                        elseif tmp==2, HDR.Patient.Height = HDR.Patient.Height*2.54; %unit='inches';
                                        elseif tmp==3, HDR.Patient.Height = HDR.Patient.Height*0.1; %unit='mm';
                                        else warning('units of height not specified');
                                        end;
                                        end;
                                elseif (tag == 7)
                                	if any(field(1:3)),
                                        HDR.Patient.Weight = field(1:2)*[1;256];
                                        tmp = field(3);
                                        if tmp==1, % unit='kg';
                                        elseif tmp==2, HDR.Patient.Weight = HDR.Patient.Weight/1000; %unit='g';
                                        elseif tmp==3, HDR.Patient.Weight = HDR.Patient.Weight/2.2; %unit='pound';
                                        elseif tmp==4, HDR.Patient.Weight = HDR.Patient.Weight*0.0284; %unit='ounce';
                                        else warning('units of weight not specified');
                                        end;
                                        end;
                                elseif tag == 8,
                                        ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
                                        HDR.Patient.Sex = field;
                                elseif tag == 9,
                                        HDR.Patient.Race = field;
                                elseif tag == 10,
					if (field(1)~=0)
	                                        HDR.Patient.Medication = field;
					else
	                                        HDR.Patient.Medication.Code = field(2:3);
						HDR.Patient.Medication = field(4:end);
					end;
                                elseif tag == 11,
                                        HDR.Patient.BloodPressure.Systolic = field*[1;256];
                                elseif tag == 12,
                                        HDR.Patient.BloodPressure.Diastolic = field*[1;256];
                                elseif tag == 13,
                                        HDR.Patient.Diagnosis = char(field);
                                elseif tag == 14,
                                        ListOfRequiredTags(ListOfRequiredTags==tag)=[];
                                        HDR.SCP1.AcquiringDeviceID = char(field);
                                        HDR.VERSION = field(15)/10;
                                elseif tag == 15,
                                        ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
                                        HDR.SCP1.AnalyisingDeviceID = char(field);
                                elseif tag == 16,
                                        HDR.SCP1.AcquiringInstitution = char(field);
                                elseif tag == 17,
                                        HDR.SCP1.AnalyzingInstitution = char(field);
                                elseif tag == 18,
                                        HDR.SCP1.AcquiringDepartment = char(field);
                                elseif tag == 19,
                                        HDR.SCP1.AnalyisingDepartment = char(field);
                                elseif tag == 20,
                                        HDR.SCP1.Physician = char(field);
                                elseif tag == 21,
                                        HDR.SCP1.LatestComfirmingPhysician = char(field);
                                elseif tag == 22,
                                        HDR.SCP1.Technician = char(field);
                                elseif tag == 23,
                                        HDR.SCP1.Room = char(field);
                                elseif tag == 24,
                                        HDR.SCP1.Emergency = field;
                                elseif tag == 25,
                                        ListOfRequiredTags(ListOfRequiredTags==tag)=[];
                                        HDR.T0(1,1:3) = [field(1:2)*[1;256],field(3:4)];
                                elseif tag == 26,
                                        ListOfRequiredTags(ListOfRequiredTags==tag)=[];
                                        HDR.T0(1,4:6) = field(1:3);
                                elseif tag == 27,
                                        HDR.Filter.HighPass = field(1:2)*[1;256]/100;
                                elseif tag == 28,
                                        HDR.Filter.LowPass = field(1:2)*[1;256]/100;
                                elseif tag == 29,
                                        if (field==0)
                                                HDR.FILTER.Notch = NaN;
                                        elseif bitand(field,1)
                                                HDR.FILTER.Notch = 60; % 60Hz Notch
                                        elseif bitand(field,2)
                                                HDR.FILTER.Notch = 50; % 50Hz Notch
                                        elseif bitand(field,3)==0
                                                HDR.FILTER.Notch = -1; % Notch Off
                                        end;
                                        HDR.SCP1.Filter.BitMap = field;
                                elseif tag == 30,
                                        HDR.SCP1.FreeText = char(field);
                                elseif tag == 31,
                                        HDR.SCP1.ECGSequenceNumber = char(field);
                                elseif tag == 32,
                                        HDR.SCP1.MedicalHistoryCodes = char(field);
                                elseif tag == 33,
                                        HDR.SCP1.ElectrodeConfigurationCodes = field;
                                elseif tag == 34,
                                        ListOfRecommendedTags(ListOfRecommendedTags==tag)=[];
                                        HDR.SCP1.Timezone = field;
                                elseif tag == 35,
                                        HDR.SCP1.MedicalHistory = char(field);
                                elseif tag == 255,
                                        % section terminator
                                elseif tag >= 200,
                                	% manufacturer specific - not standardized
                                else
                                        fprintf(HDR.FILE.stderr,'Warning SCOPEN: unknown tag %i (section 1)\n',tag);
                                end;
                        end;
                        if ~isempty(ListOfRequiredTags)
                                fprintf(HDR.FILE.stderr,'Warning SCPOPEN: the following tags are required but missing in file %s\n',HDR.FileName);
                                disp(ListOfRequiredTags);
                        end;
                        if ~isempty(ListOfRecommendedTags)
                                fprintf(HDR.FILE.stderr,'Warning SCPOPEN: the following tags are recommended but missing in file %s\n',HDR.FileName);
                                disp(ListOfRecommendedTags);
                        end;

                elseif section.ID==2, 	% Huffman tables
                        HDR.SCP2.NHT = fread(fid,1,'uint16');
                        HDR.SCP2.NCT = fread(fid,1,'uint16');
			if HDR.SCP2.NHT~=19999,
				NHT = HDR.SCP2.NHT;
			else
				NHT = 0;
			end;
                        k3 = 0;
                        for k1 = 1:NHT,
                        	HT1 = zeros(HDR.SCP2.NCT,5);
                                for k2 = 1:HDR.SCP2.NCT,
                                	tmp = fread(fid,3,'uint8') ;
                                        HDR.SCP2.prefix = tmp(1);	% PrefixLength
                                        HDR.SCP2.codelength = tmp(2);	% CodeLength
                                        HDR.SCP2.TableModeSwitch = tmp(3);
                                        tmp(4) = fread(fid,1,'int16');  % BaseValue
                                        tmp(5) = fread(fid,1,'uint32'); % BaseCode
                                	k3 = k3   + 1;
				        HT (k3,:) = [tmp'];
				        HT1(k2,:) = [tmp'];
                                end;
                                HDR.SCP2.HTree{k1} = makeTree(HT1);
                                HDR.SCP2.HTs{k1} = HT1;
			end;
			if HDR.SCP2.NHT~=19999,
				HDR.SCP2.HT = HT;
			else
				tmp = size(HT19999,1);
				HDR.SCP2.HT = [ones(tmp,1),[1:tmp]',HT19999];
                                HDR.SCP2.HTree{1} = makeTree(HT19999);
                                HDR.SCP2.HTs{1} = HT19999;
			end;

                elseif section.ID==3,
                        HDR.NS = fread(fid,1,'uint8');
                        HDR.FLAG.Byte = fread(fid,1,'uint8');
                        if ~bitand(HDR.FLAG.Byte,4)
                                fprintf(HDR.FILE.stdout,'Warning SCPOPEN: not all leads simultaneously recorded - this mode is not supported.\n');
                        end;

                        HDR.FLAG.ReferenceBeat = mod(HDR.FLAG.Byte,2);
                        %HDR.NS = floor(mod(HDR.FLAG.Byte,128)/8);
                        for k = 1:HDR.NS,
                                HDR.LeadPos(k,1:2) = fread(fid,[1,2],'uint32');
                                HDR.LeadIdCode(k,1) = fread(fid,1,'uint8');
                        end;
                        HDR.N = max(HDR.LeadPos(:))-min(HDR.LeadPos(:))+1;
                        HDR.AS.SPR = HDR.LeadPos(:,2)-HDR.LeadPos(:,1)+1;
                        HDR.SPR = HDR.AS.SPR(1);
			for k = 2:HDR.NS
				HDR.SPR = lcm(HDR.SPR,HDR.AS.SPR(k));
			end;

                        HDR = leadidcodexyz(HDR);
                        for k = 1:HDR.NS,
                                if 0,
                                elseif (HDR.LeadIdCode(k)==0),
                                        HDR.Label{k} = 'unspecified lead';
                                elseif (HDR.VERSION <= 1.3) & (HDR.LeadIdCode(k) < 86),
                                %        HDR.Label{k} = H.Label(H.LeadIdCode==HDR.LeadIdCode(k));
                                elseif (HDR.VERSION <= 1.3) & (HDR.LeadIdCode(k) > 99),
                                        HDR.Label{k} = 'manufacturer specific';
                                elseif (HDR.VERSION >= 2.0) & (HDR.LeadIdCode(k) < 151),
                                %        HDR.Label{k} = H.Label(H.LeadIdCode==HDR.LeadIdCode(k));
                                elseif (HDR.VERSION >= 2.0) & (HDR.LeadIdCode(k) > 199),
                                        HDR.Label{k} = 'manufacturer specific';
                                else
                                        HDR.Label{k} = 'reserved';
                                end;
                        end;
                        HDR.Label = strvcat(HDR.Label);

                elseif section.ID==4,
                        HDR.SCP4.L = fread(fid,1,'int16');
                        HDR.SCP4.fc0 = fread(fid,1,'int16');
                        HDR.SCP4.N = fread(fid,1,'int16');
                        HDR.SCP4.type = fread(fid,[7,HDR.SCP4.N],'uint16')'*[1,0,0,0; 0,1,0,0;0,2^16,0,0; 0,0,1,0;0,0,2^16,0; 0,0,0,1;0,0,0,2^16];

                        tmp = fread(fid,[2*HDR.SCP4.N],'uint32');
                        HDR.SCP4.PA = reshape(tmp,2,HDR.SCP4.N)';
                        HDR.SCP4.pa = [0;tmp;HDR.N];

                elseif any(section.ID==[5,6]),

                        SCP = [];
                        SCP.Cal = fread(fid,1,'int16')/1e6;    % quant in nV, converted into mV
                        SCP.PhysDim = 'mV';
                        SCP.Dur = fread(fid,1,'int16');
                        SCP.SampleRate = 1e6/SCP.Dur;
                        SCP.FLAG.DIFF  = fread(fid,1,'uint8');
                        SCP.FLAG.bimodal_compression = fread(fid,1,'uint8');

                        if isnan(HDR.NS),
				HDR.ERROR.status = -1;
				HDR.ERROR.message = sprintf('Error SCPOPEN: could not read %s\n',HDR.FileName);
				fprintf(HDR.FILE.stderr,'Error SCPOPEN: could not read %s\n',HDR.FileName);
				return;
			end;

			if CHAN==0, CHAN = 1:HDR.NS; end;
                        SCP.SPR = fread(fid,HDR.NS,'uint16');
			HDR.InChanSelect = CHAN;

                        if section.ID==6,
                                HDR.HeadLen = ftell(fid);
                                HDR.FLAG.DIFF = SCP.FLAG.DIFF;
                                HDR.FLAG.bimodal_compression = SCP.FLAG.bimodal_compression;
                                HDR.data = [];
                                outlen = HDR.SPR;
			        HDR.Calib = sparse(2:HDR.NS+1, 1:HDR.NS, SCP.Cal);
                        elseif isfield(HDR,'SCP4') %% HACK: do no know whether it is correct
                        	outlen = floor(1000*HDR.SCP4.L/SCP.Dur);
                        else
                        	outlen = inf;
                        end;

                        if ~isfield(HDR,'SCP2'),
                                if any(SCP.SPR(1)~=SCP.SPR),
                                        error('SCPOPEN: SPR do not fit');
                                else
                                        S2 = fread(fid,[SCP.SPR(1)/2,HDR.NS],'int16');
                                end;
				%S2 = S2(:,CHAN);

                        elseif (HDR.SCP2.NHT==1) && (HDR.SCP2.NCT==1) && (HDR.SCP2.prefix==0),
				codelength = HDR.SCP2.HT(1,4);
                                if (codelength==16)
                                        S2 = fread(fid,[HDR.N,HDR.NS],'int16');
                                elseif (codelength==8)
                                        S2 = fread(fid,[HDR.N,HDR.NS],'int8');
                                else
                                        fprintf(HDR.FILE.stderr,'Warning SCPOPEN: codelength %i is not supported yet.',codelength);
                                        fprintf(HDR.FILE.stderr,' Contact <a.schloegl@ieee.org>\n');
                                        return;
                                end;
				%S2 = S2(:,CHAN);

                        elseif 1, HDR.SCP2.NHT~=19999;
                        	%% User specific Huffman table
                        	%% a more elegant Huffman decoder is used here %%
                                for k = 1:HDR.NS,
                                        SCP.data{k} = fread(fid,SCP.SPR(k),'uint8');
                                end;
%                                S2 = repmat(NaN,outlen,length(HDR.InChanSelect));
                                clear S2;
                                sz = inf;
                                for k3 = 1:length(HDR.InChanSelect), k = HDR.InChanSelect(k3); %HDR.NS,
					outdata{k3} = DecodeHuffman(HDR.SCP2.HTree,HDR.SCP2.HTs,SCP.data{k},outlen);
					sz = min(sz,length(outdata{k3}));
				end;

                                for k3 = 1:length(HDR.InChanSelect), k = HDR.InChanSelect(k3); %HDR.NS,
                                	S2(:,k) = outdata{k3}(1:sz);
                                end;
				accu=0;

                        elseif HDR.SCP2.NHT==19999,
                                HuffTab = DHT;
                                for k = 1:HDR.NS,
                                        SCP.data{k} = fread(fid,SCP.SPR(k),'uint8');
                                end;
                                %for k = 1:HDR.NS,
                                for k3 = 1:length(HDR.InChanSelect), k = HDR.InChanSelect(k3); %HDR.NS,
                                %for k = CHAN(:)',
                                        s2 = SCP.data{k};
                                        s2 = [s2; repmat(0,ceil(max(HDR.SCP2.HT(:,4))/8),1)];
					k1 = 0;
					l2 = 0;
					accu = 0;
					c  = 0;
					x  = [];
					HT = HDR.SCP2.HT(find(HDR.SCP2.HT(:,1)==1),3:7);
					while (l2 < HDR.LeadPos(k,2)),
						while ((c < max(HT(:,2))) & (k1<length(s2)-1));
							k1 = k1 + 1;
							dd = s2(k1);
							accu = accu + ACC(dd+1)*(2^c);
							c = c + 8;

							if 0, %for k2 = 1:8,
								accu = accu + (dd>127)*(2^c);
								dd = mod(dd*2,256);
								c = c + 1;
							end;
						end;

                                                ixx = 1;
                                                %acc = mod(accu,2^32);   % bitand returns NaN if accu >= 2^32
						acc = accu - 2^32*fix(accu*(2^(-32)));   % bitand returns NaN if accu >= 2^32
						while (bitand(acc,2^HT(ixx,1)-1) ~= HT(ixx,5)),
							ixx = ixx + 1;
						end;

                                                dd = HT(ixx,2) - HT(ixx,1);
						if HT(ixx,3)==0,
							HT = HDR.SCP2.HT(find(HDR.SCP2.HT(:,1)==HT(ixx,5)),3:7);
							fprintf(HDR.FILE.stderr,'Warning SCPOPEN: Switching Huffman Tables is not tested yet.\n');
						elseif (dd==0),
							l2 = l2 + 1;
							x(l2) = HT(ixx,4);
						else %if (HT(ixx,3)>0),
							l2 = l2 + 1;
							%acc2  = fix(accu*(2^(-HT(ixx,1))));
							%tmp = mod(fix(accu*(2^(-HT(ixx,1)))),2^dd);

                                                        tmp = fix(accu*(2^(-HT(ixx,1))));       % bitshift(accu,-HT(ixx,1))
                                                        tmp = tmp - (2^dd)*fix(tmp*(2^(-dd)));  % bitand(...,2^dd)

                                                        %tmp = bitand(accu,(2^dd-1)*(2^HT(ixx,1)))*(2^-HT(ixx,1));
                                                        % reverse bit-pattern
                                                        if dd==8,
                                                                tmp = ACC(tmp+1);
                                                        else
                                                                tmp = dec2bin(tmp);
                                                                tmp = [char(repmat('0',1,dd-length(tmp))),tmp];
                                                                tmp = bin2dec(tmp(length(tmp):-1:1));
                                                        end
                                                        x(l2) = tmp-(tmp>=(2^(dd-1)))*(2^dd);
						end;
						accu = fix(accu*2^(-HT(ixx,2)));
						c = c - HT(ixx,2);
					end;
					x = x(:);
                                        if k3==1,
                                                S2=x(:,ones(1,k));
                                        elseif size(x,1)==size(S2,1),
                                                S2(:,k) = x;
					else
	                                        fprintf(HDR.FILE.stderr,'Error SCPOPEN: Huffman decoding failed (%i) \n',size(x,1));
	    					HDR.data = S2;
						return;
                                        end;
				end;

                        elseif (HDR.SCP2.NHT==19999), % alternative decoding algorithm.
                                warning('this branch is experimental - it might be broken')
                                HuffTab = DHT;
                                for k = 1:HDR.NS,
                                        SCP.data{k} = fread(fid,SCP.SPR(k),'uint8');
                                end;
                                %for k = 1:HDR.NS,
                                for k3 = 1:length(HDR.InChanSelect), k = HDR.InChanSelect(k3); %HDR.NS,
                                %for k = CHAN(:)',
				        tmp = SCP.data{k};
                                        accu = [tmp(4)+256*tmp(3)+65536*tmp(2)+2^24*tmp(1)];
                                        %accu = bitshift(accu,HDR.SCP2.prefix,32);
                                        c  = 0; %HDR.SCP2.prefix;
                                        l  = 4;
                                        l2 = 0;
                                        clear x;
                                        Ntmp = length(tmp);
                                        tmp = [tmp; zeros(4,1)];
                                        while c <= 32, %1:HDR.SPR(k),
                                                ixx = 1;
                                                while (bitand(accu,mask(ixx)) ~= PREFIX(ixx)),
                                                        ixx = ixx + 1;
                                                end;

                                                if ixx < 18,
                                                        c = c + prefix(ixx);
                                                        %accu  = bitshift(accu, prefix(ixx),32);
                                                        accu  = mod(accu.*(2^prefix(ixx)),2^32);
                                                        l2    = l2 + 1;
                                                        x(l2) = HuffTab(ixx,1);

                                                elseif ixx == 18,
                                                        c = c + prefix(ixx) + 8;
                                                        %accu = bitshift(accu, prefix(ixx),32);
                                                        accu  = mod(accu.*(2^prefix(ixx)),2^32);
                                                        l2    = l2 + 1;

                                                        acc1  = mod(floor(accu*2^(-24)),256);
                                                        %accu = bitshift(accu, 8, 32);
                                                        accu  = mod(accu*256, 2^32);

                                                        x(l2) = acc1-(acc1>=2^7)*2^8;
                                                        acc2  = 0;
                                                        for kk = 1:8,
                                                                acc2 = acc2*2 + mod(acc1,2);
                                                                acc1 = floor(acc1/2);
                                                        end;

                                                elseif ixx == 19,
                                                        c = c + prefix(ixx);
                                                        %accu = bitshift(accu, prefix(ixx),32);
                                                        accu  = mod(accu.*(2^prefix(ixx)),2^32);
                                                        l2    = l2 + 1;
                                                        while (c > 7) & (l < Ntmp),
                                                                l = l+1;
                                                                c = c-8;
                                                                accu = accu + tmp(l)*2^c;
                                                        end;

                                                        acc1 = mod(floor(accu*2^(-16)),2^16);
                                                        %accu = bitshift(accu, 16, 32);
                                                        accu = mod(accu.*(2^16), 2^32);

                                                        x(l2) = acc1-(acc1>=2^15)*2^16;
                                                        acc2 = 0;
                                                        for kk=1:16,
                                                                acc2 = acc2*2+mod(acc1,2);
                                                                acc1 = floor(acc1/2);
                                                        end;
                                                        %x(l2) = acc2;
                                                        c = c + 16;
                                                end;

                                                while (c > 7) & (l < Ntmp),
                                                        l = l+1;
                                                        c = c-8;
                                                        accu = accu + tmp(l)*(2^c);
                                                end;
                                        end;

                                        x = x(1:end-1)';
                                        if k3==1,
                                                S2=x(:,ones(1,k));
                                        elseif size(x,1)==size(S2,1),
                                                S2(:,k) = x;
                                        elseif 1,
	                                        fprintf(HDR.FILE.stderr,'Error SCPOPEN: length=%i of channel %i different to length=%i of channel 1 \n',size(x,1),k,size(S2,1));
                                                return;
                                        else
	                                        fprintf(HDR.FILE.stderr,'Error SCPOPEN: Huffman decoding failed (%i) \n',size(x,1));
	    					HDR.data=S2;
						return;
                                        end;
                                end;

                        elseif HDR.SCP2.NHT~=19999,
                        	%% OBSOLETE %%
                                fprintf(HDR.FILE.stderr,'Error SOPEN SCP-ECG: user specified Huffman Table not supported\n');
                                HDR.SCP = SCP;
                                return;

                        else
                                HDR.SCP2,
                        end;

                        % Decoding of Difference encoding
                        if SCP.FLAG.DIFF==2,
                                for k1 = 3:size(S2,1);
                                        S2(k1,:) = S2(k1,:) + [2,-1] * S2(k1-(1:2),:);
                                end;
                        elseif SCP.FLAG.DIFF==1,
                                S2 = cumsum(S2);
                        end;

                        if section.ID==5,
                                HDR.SCP5 = SCP;
                                HDR.SCP5.data = S2;
                                HDR.SampleRate = SCP.SampleRate;

                        elseif section.ID==6,
                                HDR.SCP6 = SCP;
                                HDR.SampleRate = SCP.SampleRate;
                                HDR.PhysDim = repmat({HDR.SCP6.PhysDim},HDR.NS,1);
                                HDR.data = S2;

                                if HDR.FLAG.bimodal_compression,
                                	%% FIXME: THIS IS A HACK - DO NOT KNOW WHETHER IT IS CORRECT.
                                %	HDR.FLAG.bimodal_compression = isfield(HDR,'SCP5') & isfield(HDR,'SCP4');
                                	HDR.FLAG.bimodal_compression = isfield(HDR,'SCP4');
				end;
                                if HDR.FLAG.bimodal_compression,
					if isfield(HDR,'SCP5')
	                                        F = HDR.SCP5.SampleRate/HDR.SCP6.SampleRate;
        	                                HDR.SampleRate = HDR.SCP5.SampleRate;
                	                        HDR.FLAG.F = F;
					else
                	                        HDR.FLAG.F = 1;
					end;

                                        tmp=[HDR.SCP4.PA(:,1);HDR.LeadPos(1,2)]-[1;HDR.SCP4.PA(:,2)+1];
                                        if ~all(tmp==floor(tmp))
                                                tmp,
                                        end;
                                        t  = (1:HDR.N) / HDR.SampleRate;
                                        S1 = zeros(HDR.N, HDR.NS);

                                        p  = 1;
                                        k2 = 1;
                                        pa = [HDR.SCP4.PA;NaN,NaN];
                                        flag = 1;
                                        %% FIXME: accu undefined ##
					accu = 0;
                                        for k1 = 1:HDR.N,
                                                if k1 == pa(p,2)+1,
                                                        flag = 1;
                                                        p    = p+1;
                                                        accu = S2(k2,:);
                                                elseif k1 == pa(p,1),
                                                        flag = 0;
                                                        k2 = ceil(k2);
                                                end;

                                                if flag,
                                                        S1(k1,:) = ((F-1)*accu + S2(fix(k2),:)) / F;
                                                        k2 = k2 + 1/F;
                                                else
                                                        S1(k1,:) = S2(k2,:);
                                                        k2 = k2 + 1;
                                                end;
                                        end;

                                        HDR.SCP.S2 = S2;
                                        HDR.SCP.S1 = S1;
                                        S2 = S1;
                                end;

                                if HDR.FLAG.ReferenceBeat & ~isfield(HDR,'SCP5')
	                                fprintf(HDR.FILE.stderr,'Warning SOPEN SCP-ECG: Flag ReferenceBeat set, but no section 5 (containing the reference beat) is available\n');
                                elseif HDR.FLAG.ReferenceBeat,

                                	tmp_data = HDR.SCP5.data*(HDR.SCP5.Cal/HDR.SCP6.Cal);
                                        for k = find(~HDR.SCP4.type(:,1)'),
                                                t1 = (HDR.SCP4.type(k,2):HDR.SCP4.type(k,4));
                                                t0 = t1 - HDR.SCP4.type(k,3) + HDR.SCP4.fc0;
                                                S2(t1,:) = S2(t1,:) + tmp_data(t0,:);
                                        end;
                                end;
	                        HDR.data  = S2;
                        end;

                elseif section.ID==7,
                        HDR.SCP7.byte1   = fread(fid,1,'uint8');
                        HDR.SCP7.Nspikes = fread(fid,1,'uint8');
                        HDR.SCP7.meanPPI = fread(fid,1,'uint16');
                        HDR.SCP7.avePPI  = fread(fid,1,'uint16');

                        for k=1:HDR.SCP7.byte1,
                                HDR.SCP7.RefBeat{k} = fread(fid,16,'uint8');
                                %HDR.SCP7.RefBeat1 = fread(fid,16,'uint8');
                        end;

                        for k=1:HDR.SCP7.Nspikes,
                                tmp = fread(fid,16,'uint16');
                                tmp(1,2) = fread(fid,16,'int16');
                                tmp(1,3) = fread(fid,16,'uint16');
                                tmp(1,4) = fread(fid,16,'int16');
                                HDR.SCP7.ST(k,:) = tmp;
                        end;
                        for k=1:HDR.SCP7.Nspikes,
                                tmp = fread(fid,6,'uint8');
                                HDR.SCP7.ST2(k,:) = tmp;
                        end;
                        HDR.SCP7.Nqrs = fread(fid,1,'uint16');
                        HDR.SCP7.beattype = fread(fid,HDR.SCP7.Nqrs,'uint8');

                        HDR.SCP7.VentricularRate = fread(fid,1,'uint16');
                        HDR.SCP7.AterialRate = fread(fid,1,'uint16');
                        HDR.SCP7.QTcorrected = fread(fid,1,'uint16');
                        HDR.SCP7.TypeHRcorr = fread(fid,1,'uint8');

                        len = fread(fid,1,'uint16');
                        tag = 255*(len==0);
                        k1 = 0;
                        while tag~=255,
                                tag = fread(fid,1,'uchar');
                                len = fread(fid,1,'uint16');
                                field = fread(fid,[1,len],'uchar');

                                if tag == 0,
                                        HDR.Patient.LastName = char(field);
                                elseif tag == 1,

                                end;
                        end;
                        HDR.SCP7.P_onset = fread(fid,1,'uint16');
                        HDR.SCP7.P_offset = fread(fid,1,'uint16');
                        HDR.SCP7.QRS_onset = fread(fid,1,'uint16');
                        HDR.SCP7.QRS_offset = fread(fid,1,'uint16');
                        HDR.SCP7.T_offset = fread(fid,1,'uint16');
                        HDR.SCP7.P_axis = fread(fid,1,'uint16');
                        HDR.SCP7.QRS_axis = fread(fid,1,'uint16');
                        HDR.SCP7.T_axis = fread(fid,1,'uint16');

                elseif section.ID==8,
                        tmp = fread(fid,9,'uint8');
                        HDR.SCP8.Report = tmp(1);
                        HDR.SCP8.Time = [[1,256]*tmp(2:3),tmp(4:8)'];
                        HDR.SCP8.N = tmp(9);
                        for k = 1:HDR.SCP8.N,
                                ix  = fread(fid,1,'uint8');
                                len = fread(fid,1,'uint16');
                                tmp = fread(fid,[1,len],'uchar');
                                HDR.SCP8.Statement{k,1} = char(tmp);
                        end

                %elseif section.ID==9,
                %        HDR.SCP9.byte1 = fread(fid,1,'uint8');

                elseif section.ID==10,
                        tmp = fread(fid,2,'uint16');
                        HDR.SCP10.NumberOfLeads = tmp(1);
                        HDR.SCP10.ManufacturerCode = tmp(2);
                        for k = []; 1:HDR.SCP10.NumberOfLeads,
                                tmp = fread(fid,2,'uint16')
                                LeadId = tmp(1);
                                LeadLen = tmp(2);
                                tmp = fread(fid,LeadLen/2,'uint16');
                                HDR.SCP10.LeadId(k)=LeadId;
                                HDR.SCP10.LeadLen(k)=LeadLen;
                                HDR.SCP10.Measurements{k}=tmp;
                        end;

                elseif section.ID==11,
                        bytes = fread(fid,11,'uint8');
                        HDR.SCP11.T0 = [bytes(2)*256+bytes(3), bytes(4:8)];
                        HDR.SCP11.Confirmed = bytes(1);
                        HDR.SCP11.NumberOfStatements = bytes(9);
                        for k = 1:HDR.SCP11.NumberOfStatements,
                                SeqNo = fread(fid,1,'uint8');
                                len11 = fread(fid,1,'uint16');
                                typeID = fread(fid,1,'uint8');
                                Statement = fread(fid,len11-1,'uint8');
                                HDR.SCP11.Statement.SeqNo(k) = SeqNo;
                                HDR.SCP11.Statement.len11(k) = len11;
                                HDR.SCP11.Statement.typeID(k) = typeID;
                                HDR.SCP11.Statement.Statement{k} = Statement;
                        end;
                end;

		if ~section.Length,
			HDR.ERROR.status  = -1;
			HDR.ERROR.message = 'Error SCPOPEN: \n';
			return;
		end;
        end;

        HDR.SPR  = size(HDR.data,1);
        HDR.NRec = 1;
        HDR.AS.endpos = HDR.SPR;

        HDR.FILE.OPEN = 0;
        HDR.FILE.POS  = 0;
        HDR.TYPE = 'native';
        fclose(HDR.FILE.FID);

else    % writing SCP file

	NSections = 12;
        SectIdHdr = zeros(1,16);
        VERSION = round(HDR.VERSION*10);
        if ~any(VERSION==[10,13,20])
                fprintf(HDR.FILE.stderr,'Warning SCPOPEN(WRITE): unknown Version number %4.2f\n',HDR.VERSION);
              	VERSION = 20;
        end;
	SectIdHdr(9:10) = VERSION; % Section and Protocol version number

        POS = 6; B = zeros(1,POS);
        for K = 0:NSections-1;
                b = [];
                if K==0,
                        % SECTION 0
                        b = [SectIdHdr(1:10),'SCPECG', zeros(1,NSections*10)];
	                b(16+(7:10)) = s4b(POS+1);

                elseif K==1,
                        % SECTION 1
                        b = SectIdHdr;
                        % tag(1),len(1:2),field(1:len)
                        if isfield(HDR.Patient,'Name'),
				b = [b, 0, s2b(length(HDR.Patient.Name)), HDR.Patient.Name];
			end;
                        if isfield(HDR.Patient,'Id'),
				b = [b, 2, s2b(length(HDR.Patient.Id)), HDR.Patient.Id];
			end;
                        if isfield(HDR.Patient,'Age'),
				%b = [b, 4, s2b(3), s2b(HDR.Patient.Age),1];  %% use birthday instead
			end;
                        if isfield(HDR.Patient,'Birthday'),
	                        b = [b, 5, s2b(4), s2b(HDR.Patient.Birthday(1)),HDR.Patient.Birthday(2:3)];
			end;
                        if isfield(HDR.Patient,'Height'),
                        if ~isnan(HDR.Patient.Height),
                        	b = [b, 6, s2b(3), s2b(HDR.Patient.Height),1];
                        end;
                        end;
                        if isfield(HDR.Patient,'Weight'),
                        if ~isnan(HDR.Patient.Weight),
                        	b = [b, 7, s2b(3), s2b(HDR.Patient.Weight),1];
                        end;
                        end;
	                if isfield(HDR.Patient,'Sex'),
	                if ~isempty(HDR.Patient.Sex)
                        	sex = HDR.Patient.Sex;
                        	if 0,
                        	elseif isnumeric(sex),           sex = sex(1);
                        	elseif strncmpi(sex,'male',1);   sex = 1;
                        	elseif strncmpi(sex,'female',1); sex = 2;
                        	else sex = 9; % unspecified
                        	end;
	                      	b = [b, 8, s2b(1), sex];
	                end;
                        end;
                        if isfield(HDR.Patient,'Race'),
 	                       	b = [b, 9, s2b(1), HDR.Patient.Race(1)];
                        end;
                        if isfield(HDR.Patient,'BloodPressure'),
                        	b = [b, 11, s2b(2), s2b(HDR.Patient.BloodPressure.Systolic)];
                        	b = [b, 12, s2b(2), s2b(HDR.Patient.BloodPressure.Diastolic)];
                        end;

                        %% Tag 14
                        tag14.AnalyzingProgramRevisionNumber = ['',char(0)];
                        tag14.SerialNumberAcqDevice = ['',char(0)];
                        tag14.AcqDeviceSystemSoftware = ['',char(0)];
                        tag14.SCPImplementationSoftware = ['BioSig4OctMat v 1.76+',char(0)];
                        tag14.ManufactureAcqDevice = ['',char(0)];
                        t14 = [zeros(1,35), length(tag14.AnalyzingProgramRevisionNumber),tag14.AnalyzingProgramRevisionNumber,tag14.SerialNumberAcqDevice,tag14.AcqDeviceSystemSoftware,tag14.SCPImplementationSoftware,tag14.ManufactureAcqDevice];
                        t14(8)  = 255;  % Manufacturer
                        %%% ### FIXME ###  t14(9:14) = % cardiograph model
                        t14(15) = VERSION;        % Version
                        t14(16) = hex2dec('A0');  % Demographics and ECG rhythm data" (if we had also the reference beats we should change it in 0xC0).
                        t14(18) = hex2dec('D0');  % Capabilities of the ECG Device: 0xD0 (acquire, print and store).
                        b   = [b, 14, s2b(length(t14)), t14];

                        b = [b, 25, s2b(4), s2b(HDR.T0(1)),HDR.T0(2:3)];
                        b = [b, 26, s2b(3), HDR.T0(4:6)];
                        if ~any(isnan(HDR.Filter.HighPass))
	                        b = [b, 27, s2b(2), s2b(round(HDR.Filter.HighPass(1)*100))];
                        end;
                        if ~any(isnan(HDR.Filter.LowPass))
	                	b = [b, 28, s2b(2), s2b(round(HDR.Filter.LowPass(1)))];
			end;
                        b = [b, 255, 0, 0];	% terminator
			b = b + (b<0)*256;

                elseif K==3,
                        % SECTION 3
                        b = [SectIdHdr,HDR.NS,4+HDR.NS*8];
                        if ~isfield(HDR,'LeadIdCode'), HDR.LeadIdCode = zeros(1,HDR.NS); end;
                        if (numel(HDR.LeadIdCode)~=HDR.NS); warning('HDR.LeadIdCode does not have HDR.NS elements'); end;
                        if any(HDR.LeadIdCode>255), warning('invalid LeadIdCode'); end;
                        for k = 1:HDR.NS,
                                b = [b, s4b(1), s4b(HDR.SPR*HDR.NRec), mod(HDR.LeadIdCode(k),256)];
                        end;

                elseif K==6,
                        % SECTION 6
                        Cal = full(HDR.Calib(2:end,:)); Cal = Cal - diag(diag(Cal));
                        if any(Cal(:))
                                fprintf(HDR.FILE.stderr,'Calibration is not a diagonal matrix.\n\tThis can result in incorrect scalings.\n');
                        end;
                        Cal = full(diag(HDR.Calib(2:end,:)));
                        if any(Cal~=Cal(1)),
                                fprintf(HDR.FILE.stderr,'scaling information is not equal for all channels; \n\tThis is not supported by SCP and can result in incorrect scalings.\n');
                        end;
                        [tmp,scale1] = physicalunits(HDR.PhysDim{1});
                        [tmp,scale2] = physicalunits('nV');
                        b = [SectIdHdr, s2b(round(Cal(1)*scale1/scale2)), s2b(round(1e6/HDR.SampleRate)), 0, 0];
                        for k = 1:HDR.NS,
                                b = [b, s2b(HDR.SPR*HDR.NRec*2)];
                        end;
                        data = HDR.data(:);
                        data = data + (data<0)*2^16;
                        tmp  = s2b(round(data))';
                        b = [b,tmp(:)'];
                else
                        b = [];
                end;
                if (length(b)>0)
                        if mod(length(b),2), % align to multiple of 2-byte blocks
                                b = [b,0];
                        end;
                        if (length(b)<16), fprintf(HDR.FILE.stderr,'section header %i less then 16 bytes %i', K,length(b)); end;
                        b(3:4) = s2b(K);
                        b(5:8) = s4b(length(b));
                        %b(1:2)= s4b(crc);
                        b(1:2) = s2b(crc16eval(b(3:end)));
                        % section 0: startpos in pointer field
	                B(22+K*10+(7:10)) = s4b(POS+1);
                end;
                B = [B(1:POS),b];
                % section 0 pointer field
                B(22+K*10+(1:2)) = s2b(K);
                B(22+K*10+(3:6)) = s4b(length(b)); % length
                POS = POS + length(b);
        end
        B(3:6) = s4b(POS);      % length of file
        B(7:8) = s2b(crc16eval(B(9:22+NSections*10)));  % CRC of Section 0

        B(1:2) = s2b(crc16eval(B(3:end)));

        %        fwrite(fid,crc,'int16');
        count = fwrite(fid,B,'uchar');
        fclose(fid);
end
end  %% scpopen

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Auxillary functions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function b2 = s2b(i)
	% converts 16bit into 2 bytes
	b2 = [bitand(i,255),bitand(bitshift(i,-8),255)];
	return;
end	%%%%% s2b %%%%%

function b4 = s4b(i)
	% converts 32 bit into 4 bytes
	b4 = [s2b(bitand(i,2^16-1)),s2b(bitand(bitshift(i,-16),2^16-1)) ];
	return;
end	%%%%% s4b %%%%%

function T = makeSubTree(T,bc,len,val)
	if (len==0)
		T.idxTable = val;
		return
	end;
if 1,
	b = bitand(bc,1)+1;
	if ~isfield(T,'branch') T.branch = {[],[]}; end;
	T.branch{b} = makeSubTree(T.branch{b},bitshift(bc,-1),len-1,val);
else
	if bitand(bc,1)
		if ~isfield(T,'node1'),T.node1 = []; end;
		T.node1 = makeSubTree(T.node1,bitshift(bc,-1),len-1,val);
	else
		if ~isfield(T,'node0'),T.node0 = []; end;
		T.node0 = makeSubTree(T.node0,bitshift(bc,-1),len-1,val);
	end;
end,
	return;
end	%%%%% makeSubTree %%%%%

function T = makeTree(HT)
	T = [];
	for k1 = 1:size(HT,1)
		for k2 = 1:HT(k1,1) % CodeLength
			T = makeSubTree(T,HT(k1,5),HT(k1,1),k1);
		end;
	end;
%save matlab T,pause
	return;
end	%%%%% makeTree %%%%%

function outdata = DecodeHuffman(HTrees,HTs,indata,outlen)
	ActualTable = 1;
	k1 = 1; r=0;
	k2 = 0;

	if ((outlen>0) && isfinite(outlen))
		outdata = repmat(NaN,outlen,1);
	else
		outdata = [0;0];  %% make it a column vector
	end;
	Node    = HTrees{ActualTable};
	while ((k1*8+r <= 8*length(indata)) && (k2<outlen))
		if ~isfield(Node,'idxTable')
			r = r+1; if (r>8), k1=k1+1; r=1; end;
if 1,
			b = bitand(bitshift(indata(k1),r-8),1)+1;
			if ~isempty(Node.branch{b})
				Node = Node.branch{b};
			else
				fprintf(2,'Warning SCPOPEN: empty node in Huffman table\n');
			end;
else
			if bitand(bitshift(indata(k1),r-8),1)
				if isfield(Node,'node1')
					Node = Node.node1;
				else
					fprintf(2,'Warning SCPOPEN: empty node in Huffman table\n');
				end;
			else
				if isfield(Node,'node0')
					Node = Node.node0;
				else
					fprintf(2,'Warning SCPOPEN: empty node in Huffman table\n');
				end;
			end;
end;
		end;

		if isfield(Node,'idxTable')
			TableEntry = HTs{ActualTable}(Node.idxTable,:);
			dlen = TableEntry(2)-TableEntry(1);
			if (~TableEntry(3))
				ActualTable = TableEntry(4);
			elseif (dlen~=0)
				acc = 0;
				for k3 = 1:dlen,
					r = r+1; if (r>8), k1=k1+1; r=1; end;
					acc = 2*acc + bitand(bitshift(indata(k1),r-8), 1);
				end;
				if (acc>=bitshift(1,dlen-1))
					acc = acc - bitshift(1,dlen);
				end;
				k2 = k2+1;
				outdata(k2) = acc;
			else
				k2 = k2+1;
				outdata(k2)=TableEntry(4);
			end;
			Node = HTrees{ActualTable};
		end;
	end;
	return;
end %%%%%%%% DecodeHuffman %%%%%%%%%

function crc16 = crc16eval(D)
% CRC16EVAL cyclic redundancy check with the polynomiaL x^16+x^12+x^5+1
% i.e. CRC-CCITT http://en.wikipedia.org/wiki/Crc16

	D = uint16(D);

	crchi = 255;
	crclo = 255;

	t = '00102030405060708191a1b1c1d1e1f112023222524272629383b3a3d3c3f3e32434041464744454a5b58595e5f5c5d53626160676665646b7a79787f7e7d7c74858687808182838c9d9e9f98999a9b95a4a7a6a1a0a3a2adbcbfbeb9b8bbbab6c7c4c5c2c3c0c1cedfdcdddadbd8d9d7e6e5e4e3e2e1e0effefdfcfbfaf9f8f9181b1a1d1c1f1e110003020504070608393a3b3c3d3e3f30212223242526272b5a59585f5e5d5c53424140474645444a7b78797e7f7c7d72636061666764656d9c9f9e99989b9a95848786818083828cbdbebfb8b9babbb4a5a6a7a0a1a2a3afdedddcdbdad9d8d7c6c5c4c3c2c1c0cefffcfdfafbf8f9f6e7e4e5e2e3e0e1e';
	crc16htab = hex2dec(reshape(t,2,length(t)/2)');

	t = '0021426384a5c6e708294a6b8cadceef31107352b594f7d639187b5abd9cffde62432001e6c7a4856a4b2809eecfac8d53721130d7f695b45b7a1938dffe9dbcc4e586a740610223cced8eaf48690a2bf5d4b79671503312fddcbf9e79583b1aa687e4c522036041ae8feccd2a0b684997b6d5f4133251709fbeddfc1b3a597888a9caeb0c2d4e6f80a1c2e304254667b998fbda3d1c7f5eb190f3d235147756eacba8896e4f2c0de2c3a08166472405dbfa99b85f7e1d3cd3f291b0577615344c6d0e2fc8e98aab44650627c0e182a37d5c3f1ef9d8bb9a75543716f1d0b3922e0f6c4daa8be8c926076445a283e0c11f3e5d7c9bbad9f81736557493b2d1f0';
	crc16ltab = hex2dec(reshape(t,2,length(t)/2)');

	for k = 1:length(D),
		ix = double(bitxor(crchi,D(k)))+1;
		crchi = bitxor(crclo,crc16htab(ix));
		crclo = crc16ltab(ix);
	end;
	crc16 = crchi*256+crclo;

end	%%%%% crc16eval %%%%%

author: PatternRecognition

repo: OpenBMI


text:
function S = getcovariance(mix,inds,centers)

if ~exist('inds','var'),
  inds = 1:mix.nin;
end;
if ~exist('centers','var'),
  centers = 1:mix.ncentres;
end;

inds = inds(:)';
ninds = length(inds);
centers = centers(:)';
ncenters = length(centers);
S = zeros(ninds,ninds,ncenters);

switch mix.covar_type,

 case 'spherical',
  for i = 1:ncenters,
    j = centers(i);
    S(:,:,i) = eye(ninds)*mix.covars(j);
  end;

 case 'diag',
  for i = 1:ncenters,
    j = centers(i);
    S(:,:,i) = diag(mix.covars(j,:));
  end;

 case {'full','sparse'},
  S = mix.covars(inds,inds,centers);

 case 'ppca',
  for i = 1:ncenters,
    j = centers(i);
    covars = mix.covars(j) * eye(mix.nin) + ...
	     mix.U(:, :, j)*(diag(mix.lambda(j, :))-diag(mix.covars(j)))* ...
	     (mix.U(:, :, j)');
    S(:,:,i) = covars(inds,inds);
  end;

 case 'zero',
  S(:) = 0;

 otherwise,
  error(['Unknown covar_type ',mix.covar_type]);
end;

author: kristinbranson

repo: JAABA


text:
function seats=apport(votes,num,method,minseats)
%APPORT Political Apportionment
%   SEATS = APPORT(VOTES,NUM,METHOD) takes the array of VOTES for each
%   member and allocates seats according to METHOD until NUM seats have
%   been allocated.
%
%   APPORT(...,MINSEATS) pre-allocates MINSEATS to each member first,
%   guaranteeing that each member has at least MINSEATS number of SEATS.
%   MINSEATS can either be a vector, with length equal to the number of
%   members, or a scalar which is the same for all members. (Default is 0)
%
%   If either NUM or METHOD is multivalued it will enumerate all
%   combinations within those parameters. An error will be given if both
%   are multivalued.
%
%   Below is a list of available methods. METHOD can either be written out
%   in full ('Greatest Divisor') or in abbreviated form ('GD'). Also listed
%   are common names for the same method.
%
%   Methods available:
%   ------------------
%   Greatest Divisor   : d'Hondt; Jefferson; Highest Average
%   Smallest Divisor   : Adams
%   Equal Proportions  : Huntington; Hill; Geometric Mean
%   Harmonic Mean      : Dean
%   Major Fractions    : Sainte-Lague; Webster; Arithmetic Mean; Odd Numbers
%   Largest Remainders : Vinton; Hamilton; Hare
%
%   Abbreviated form:
%   -----------------
%   ALL : All six unique methods
%   A   : Adams
%   AM  : Arithmetic Mean
%   D   : Dean
%   DH  : d'Hont
%   EP  : Equal Proportions
%   GD  : Greatest Divisor
%   GM  : Geometric Mean
%   HH  : Huntington OR Hill (equivalent)
%   HA  : Highest Average
%   HAM : Hamilton OR Hare (equivalent)
%   HM  : Harmonic Mean
%   J   : Jefferson
%   LR  : Largest Remainders
%   MF  : Major Fractions
%   ON  : Odd Numbers
%   SD  : Smallest Divisor
%   SL  : Sainte-Lague
%   V   : Vinton
%   W   : Webster
%
%
%   EXAMPLE 1
%   Allocate 80 seats among four parties with votes of [100 80 30 20], with
%   minimum of 3 seats each. Compare all methods.
%      apport([100 80 30 20],80,'All',3)
%
%   EXAMPLE 2
%   Compare method of Greatest Divisor with the method of Smallest Divisor
%   of allocating 15 seats among four parties with votes of [104 62; 34 55]
%       apport([104 62; 34 55],15,{'GD','SD'})
%
%   EXAMPLE 3
%   'Alabama Paradox'
%   "increasing the total number of items would decrease one of the shares"
%      apport([2 6 6],10:13,'Hamilton')
%   Notice first party decreased representation from 2 to 1, even with
%   an increase in total number of shares
%
%   EXAMPLE 4
%   'New states paradox'
%   "it is possible for an existing state to get more representatives than
%   if the new state were not added"
%      apport([27 39 68  ],23,'Largest Remainders')
%      apport([27 39 68 6],23,'Largest Remainders')
%   Notice first state increased representation from 4 to 5, at the
%   expense of one of the other two original states
%

%   Mike Sheppard
%   Last Modified: 1-Feb-2012

%%ERROR CHECKING
if nargin<3
    error('APPORT:TooFewInputs','Too few inputs');
end
if nargin==3
    minseats=0;
end
sz=size(votes); votes=votes(:); minseats=minseats(:);
if isscalar(minseats)
    minseats=minseats.*ones(size(votes)); %make into vector
end
if numel(minseats)~=numel(votes)
    error('APPORT:Number','MINSEATS must be either a scalar or vector of length equal to VOTES');
end
if any(minseats<0)||any(rem(minseats,1))||(any(~isfinite(minseats)))
    error('APPORT:Number','MINSEATS must be a non-negative integers');
end
if (num<0)|(rem(num,1))|(any(~isfinite(num)))
    error('APPORT:Number','NUM must be a positive integer');
end
if (~isnumeric(votes))||(any(~isfinite(votes)))||(any(isnan(votes)))||(~isvector(votes))
    error('APPORT:Votes','VOTES must be a finite numeric value, either scalar or vector');
end
if sum(minseats)>num
    error('APPORT:Number','ERROR: Total minimum allocated seats exceed total seats allowed');
end
method=lower(method);
if ~any(ismember(method,{'all','greatest divisor','dhondt','jefferson',...
        'highest average','smallest divisor','adams','equal proportions',...
        'huntington','hill','geometric mean','harmonic mean','dean',...
        'major fractions','sainte-lague','webster','arithmetic mean',...
        'odd numbers','largest remainders','vinton','hamilton','hare',...
        'a','am','d','dh','ep','gd','gm','hh','ha','ham','hm','j','lr',...
        'mf','on','sd','sl','v','w'}))
    error('APPORT:Method','Method must one of the items in the list under ''help apport''');
end
if iscell(method)&&any(strcmp(method,'all'))
    error('APPORT:Method','METHOD of ALL must be given as a string with no other methods listed');
end
if ischar(method)&&strcmp(method,'all')
    methodv=lower({'Greatest Divisor','Smallest Divisor',...
        'Equal Proportions','Harmonic Mean','Major Fractions',...
        'Largest Remainders'});
elseif ~iscell(method)
    methodv={method}; %Singular method
else
    methodv=method; %Already a cell
end
if all([numel(methodv) numel(num)]>1)
    error('APPORT:Method','NUM and METHOD cannot both be multivalued');
end

nummethods=numel(methodv);
numnum=numel(num);
methodsingular=(nummethods==1);
maxcount=max([nummethods numnum]);

%Save original values
isvec=isvector(ones(sz));

for k=1:maxcount
    if methodsingular
        %methods is singular, loop through num
        seats=apport_helper(votes,num(k),methodv{1},minseats);
    else
        %num is singular, loop through method
        seats=apport_helper(votes,num(1),methodv{k},minseats);
    end
    seats=reshape(seats,sz);

    if maxcount==1
        output=seats; %singular input
    else
        %Multiple inputs
        if ~isvec
            %Multi-dimensional array
            temp(1:ndims(seats))={':'};
            output(temp{:},k)=seats;
        else
            %Vector
            if sz(1)==1
                output(k,:)=seats;
            else
                output(:,k)=seats;
            end
        end
    end

end

%Rename for output
seats=output;

end

function seats=apport_helper(votes,num,method,minseats)
%Subtract minseats from original number
num=num-sum(minseats);

if ismember(method,{'largest remainders','vinton','hamilton','hare','lr','v','ham'})
    %Method of Largest Remainders is handled as specific case
    st_quota=votes./(sum(votes)./num);                                %Standard Quota
    [frac_lost,indx]=sort(rem(st_quota-floor(st_quota),1),'descend'); %Sort by remainders
    seats=floor(st_quota);                                            %Allocate the integer values
    num_missing=num-sum(seats);                                       %Determine number of seats left to allocate
    seats(indx(1:num_missing))=seats(indx(1:num_missing))+1;          %Increase by one seat those with largest remainders
    seats=seats(:);                                                   %Consistent output, column vector
else            %One of the other methods

    %   Algorithm:
    %   Party A with vote total P(A) is entitled to its m'th seat before
    %   party B with vote total P(B) is entitled to its n'th seat
    %   if and only if P(A)/f(m-1) > P(B)/f(n-1) where f(x) is the function
    %   f(x) = x (Smallest Divisor)
    %   f(x) = harmonic mean of x and x + 1 (Harmonic Mean)
    %   f(x) = geometric mean of x and x + 1 (Equal Proportions)
    %   f(x) = arithmetic mean of x and x + 1 (Major Fractions)
    %   f(x) = x + 1 (Greatest Divisor)

    x=[1:num];  %Seat number
    %Note: fd=f(x-0)
    switch lower(method)
        case {'smallest divisor','adams','sd','a'}
            fd=(x);
        case {'harmonic mean','dean','hm','d'}
            fd=harmmean([(x);(x)+1]);
        case {'equal proportions','huntington','hill','geometric mean','ep','hh','gm'}
            fd=geomean([(x);(x)+1]);
        case {'major fractions','sainte-lague','webster', 'arithmetic mean',...
                'odd numbers','mf','sl','w','am','on'}
            fd=mean([(x);(x)+1]);
        case {'greatest divisor','dhondt','jefferson','highest average',...
                'gd','dh','j','ha'}
            fd=(x)+1;
        otherwise
            error('APPORT:Method','Method must one of the items in the list under ''help''');
    end

    %Allocate seats according to priority values
    pv=kron(votes(:),1./fd);               %priority value matrix
    [I,J] = ind2sub(size(pv),1:numel(pv)); %I party index of votes
    [Y,indx] = sort(pv(:),'descend');      %indx is linear index
    %Take the top 'num' of priority values and determine the
    %appropriate party index for each, and group to determine total number of
    %seats for each party
    list=I(indx(1:num));
    seats=accumarray(list',ones(size(list)));

    %Special case to catch:
    %zero seats at end do not show up, add if necessary
    if length(seats)<length(votes)
        seats=[seats; zeros(length(votes)-length(seats),1)];
    end

end

%Add back minimum number of seats and reshape to original shape
seats=seats+minseats;
end

author: Sable

repo: mcbench-benchmarks


text:
function [e0,e1,theta] = det2MM(im,sigmaI,sigmaO)
% function [e0,e1,theta] = det2MM(im,sigmaI,sigmaO)
%
% Compute the eigenspectrum of the spatially averaged second moment
% matrix.
%
% INPUT
%	im	Image.
%	sigmaI	Inner scale (sigma for image derivatives).
%	sigmaO	Outer scale (sigma for spatial averaging).
%
% OUTPUT
%	e0,e1	Smaller,larger eigenvalues.
%	theta	Orientation of fist eigenvector + pi/2
%		(i.e. orientation of possible edge).
%
% David R. Martin <dmartin@eecs.berkeley.edu>
% March 2003

if ndims(im)>2, im = rgb2gray(im); end;
idiag = norm(size(im));

if nargin<2, sigmaI=2; end
if nargin<3, sigmaO=sigmaI; end % Scott Konishi says this is close
                                % to optimal.
sigmaI = max(0.5,sigmaI);
sigmaO = max(0.5,sigmaO);

% compute x and y image derivatives at inner scale
fb = cell(2,1);
fb{1} = oeFilter(sigmaI,3,pi/2,1);
fb{2} = fb{1}';
fim = fbRun(fb,im);
dx = fim{1};
dy = fim{2};

% compute smoothed squared image derivatives at outer scale
f = oeFilter(sigmaO,3);
dx2 = applyFilter(f,dx.^2);
dy2 = applyFilter(f,dy.^2);
dxy = applyFilter(f,dx.*dy);

% compute eigenvalues of the spatially averaged 2nd moment matrix
% and the orientations of the eigenvectors
k = sqrt( (dx2-dy2).^2 + 4.*dxy.^2 );
eig0 = (dx2 + dy2 - k) / 2;
eig1 = (dx2 + dy2 + k) / 2;
t0 = atan2( dx2-eig0, -dxy );
t1 = atan2( dx2-eig1, -dxy );

% order eigenvalues by their absolute value, so e0<=e1, and pick
% out the orientation corresponding to the largest eigenvalue
x = (abs(eig1) > abs(eig0));
e0 = abs(eig0.*x + eig1.*~x);
e1 = abs(eig1.*x + eig0.*~x);
theta = t1.*x + t0.*~x;
theta = mod(theta+pi/2,pi);

% check postconditions
if any(e0>e1), error('e0>e1'); end

author: Cloud-CV

repo: object-proposals


text:
% Test file for periodic eigenvalue problems.
function pass = test_eigs_periodic(pref)

% Get the preferences.
if ( nargin < 1 )
    pref = cheboppref;
end

tol = 1e2*pref.bvpTol;

%% EXAMPLE USING TRIGCOLLOC.
% Solving
%   -u'' = lambda*u,
% for x in [0 2*pi],  subject to
% periodic boundary conditions.

% Define the domain we're working on.
dom = [0, 2*pi];

% Assign the equation to a chebop N such that N(u) = lambda*u.
L = chebop(@(u) -diff(u, 2), dom);

% Assign boundary conditions to the chebop.
L.bc = 'periodic';

% Number of eigenvalue and eigenmodes to compute.
k = 7;

% Solve the eigenvalue problem with FOURIER.
[V, D] = eigs(L, k);
D = diag(D);

% Exact solution.
Dexact = [0 1 1 4 4 9 9]';
pass(1) = norm(D - Dexact, inf) < tol;
pass(2) = isequal(get(V{1}.funs{1}, 'tech'), @trigtech);

%% ANOTHER EXAMPLE USING TRIGCOLLOC.
% Solving
%   -u'' + 2*q*cos(2*x)*u = lambda*u,
% for x in [0 2*pi],  subject to
% periodic boundary conditions.

% Define the domain we're working on.
dom = [0, 2*pi];

% Assign the equation to a chebop N such that N(u) = lambda*u.
q = 2;
L = chebop(@(x, u) -diff(u, 2) + 2*q*cos(2*x).*u, dom);

% Assign boundary conditions to the chebop.
L.bc = 'periodic';

% Number of eigenvalue and eigenmodes to compute.
k = 7;

% Solve the eigenvalue problem with FOURIER.
[V, D] = eigs(L, k);
D = diag(D);

% Solutions from WolframAlpha.
Dwolfram = [ -1.513956885056520;
             -1.390676501225323;
              2.379199880488686;
              3.672232706497191;
              5.172665133358294;
              9.140627737766440;
              9.370322483621104 ];
pass(3) = norm(D - Dwolfram, inf) < tol;
pass(4) = isequal(get(V{1}.funs{1}, 'tech'), @trigtech);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function result = cone_unit_3d ( func )

%*****************************************************************************80
%
%% CONE_UNIT_3D approximates an integral inside a unit cone in 3D.
%
%  Integration Region:
%
%    X**2 + Y**2 <= 1 - Z
%    0 <= Z <= 1.
%
%  Discussion:
%
%    An 48 point degree 7 formula, Stroud CN:S2:7-1, is used.
%
%    (There is a typographical error in the S2:7-1 formula for B3.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971, page 339.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied function which
%    evaluates F(X,Y,Z), of the form
%      function value = func ( x, y, z )
%
%    Output, real RESULT, the approximate integral of the function.
%
  u = [ 0.04850054945E+00, 0.2386007376E+00, 0.5170472951E+00,  0.7958514179E+00 ];
  w1 = [ 0.1108884156E+00,  0.1434587878E+00, 0.06863388717E+00, 0.01035224075E+00 ];

  a = sqrt ( 3.0E+00 ) / 2.0E+00;
  b = sqrt ( ( 27.0E+00 - 3.0E+00 * sqrt ( 29.0E+00 ) ) / 104.0E+00 );
  c = sqrt ( ( 27.0E+00 + 3.0E+00 * sqrt ( 29.0E+00 ) ) / 104.0E+00 );
  w2(1:3) = 3.0E+00 * [ ...
    2.0E+00 / 27.0E+00, ...
    ( 551.0E+00 + 4.0E+00 * sqrt ( 29.0E+00 ) ) / 6264.0E+00, ...
    ( 551.0E+00 - 4.0E+00 * sqrt ( 29.0E+00 ) ) / 6264.0E+00 ];

  quad = 0.0E+00;

  for i = 1 : 4

    x = a * ( 1.0E+00 - u(i) );
    y = 0.0E+00;
    z = u(i);
    quad = quad + w1(i) * w2(1) * feval ( func, x, y, z );

    x = - a * ( 1.0E+00 - u(i) );
    y = 0.0E+00;
    z = u(i);
    quad = quad + w1(i) * w2(1) * feval ( func, x, y, z );

    x = 0.0E+00;
    y = a * ( 1.0E+00 - u(i) );
    z = u(i);
    quad = quad + w1(i) * w2(1) * feval ( func, x, y, z );

    x = 0.0E+00;
    y = - a * ( 1.0E+00 - u(i) );
    z = u(i);
    quad = quad + w1(i) * w2(1) * feval ( func, x, y, z );

  end

  for i = 1 : 4

    x =   b * ( 1.0E+00 - u(i) );
    y =   b * ( 1.0E+00 - u(i) );
    z =   u(i);
    quad = quad + w1(i) * w2(2) * feval ( func, x, y, z );

    x = - b * ( 1.0E+00 - u(i) );
    y =   b * ( 1.0E+00 - u(i) );
    z =   u(i);
    quad = quad + w1(i) * w2(2) * feval ( func, x, y, z );

    x = - b * ( 1.0E+00 - u(i) );
    y = - b * ( 1.0E+00 - u(i) );
    z =   u(i);
    quad = quad + w1(i) * w2(2) * feval ( func, x, y, z );

    x =   b * ( 1.0E+00 - u(i) );
    y = - b * ( 1.0E+00 - u(i) );
    z =   u(i);
    quad = quad + w1(i) * w2(2) * feval ( func, x, y, z );

    x =   c * ( 1.0E+00 - u(i) );
    y =   c * ( 1.0E+00 - u(i) );
    z =   u(i);
    quad = quad + w1(i) * w2(3) * feval ( func, x, y, z );

    x = - c * ( 1.0E+00 - u(i) );
    y =   c * ( 1.0E+00 - u(i) );
    z =   u(i);
    quad = quad + w1(i) * w2(3) * feval ( func, x, y, z );

    x = - c * ( 1.0E+00 - u(i) );
    y = - c * ( 1.0E+00 - u(i) );
    z =   u(i);
    quad = quad + w1(i) * w2(3) * feval ( func, x, y, z );

    x =   c * ( 1.0E+00 - u(i) );
    y = - c * ( 1.0E+00 - u(i) );
    z =   u(i);
    quad = quad + w1(i) * w2(3) * feval ( func, x, y, z );

  end

  r = 1.0E+00;
  h = 1.0E+00;

  volume = cone_volume_3d ( r, h );
  result = quad * volume;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [problemTypeParams, solverParams] = parseSolverParameters(problemType, varargin)
% Gets default cobra solver parameters for a problem of type problemType, unless
% overridden by cobra solver parameters provided by varagin either as parameter
% struct, or as parameter/value pairs.
%
% USAGE:
%    [problemTypeParams, solverParams] = parseSolverParameters(problemType,varargin)
%
% INPUT:
%    problemType:       The type of the problem to get parameters for
%                       ('LP','MILP','QP','MIQP','NLP','EP')
%
% OPTIONAL INPUTS:
%    varargin:          Additional parameters either as parameter struct, or as
%                       parameter/value pairs. A combination is possible, if
%                       the parameter struct is either at the beginning or the
%                       end of the optional input.
%                       All fields of the struct which are not COBRA parameters
%                       (see `getCobraSolverParamsOptionsForType`) for this
%                       problem type will be passed on to the solver in a
%                       solver specific manner.
%
% OUTPUTS:
%    problemTypeParams: The COBRA Toolbox specific parameters for this
%                       problem type given the provided parameters
%
%    solverParams:      Additional parameters provided which are not part
%                       of the COBRA parameters and are assumed to be part
%                       of direct solver input structs. For some solvers, it
%                       is essential to not include any extraneous fields that are
%                       outside the solver interface specification.

cobraSolverParameters = getCobraSolverParamsOptionsForType(problemType); % build the default Parameter Structure

% set the solver Type
eval(['global CBT_' problemType '_SOLVER;'])
eval(['defaultSolver = CBT_' problemType '_SOLVER;']);

% initialize the solver variables
solverVars = cell(numel(cobraSolverParameters),1);

% get the default variables for the correct solver.
[solverVars{:}] = getCobraSolverParams(problemType,cobraSolverParameters,struct('solver',defaultSolver));
defaultParams = [cobraSolverParameters',solverVars];

nVarargin = numel(varargin);
% parse the supplied parameters
if nVarargin > 0
    % we should have a struct at the end
    if mod(nVarargin,2) == 1
        optParamStruct = varargin{end};
        if ~isstruct(optParamStruct)
            % but it could also be at the first position, so test that as well.
            optParamStruct = varargin{1};
            varargin(1) = [];
            nVarargin = numel(varargin); %added this in case varagin{1} is the parameter structure
            if ~isstruct(optParamStruct)
                error(['Invalid Parameters supplied.\n',...
                       'Parameters have to be supplied either as parameter/Value pairs, or as struct.\n',...
                       'A combination is possible, if the last or first input argument is a struct, and all other arguments are parameter/value pairs'])
            end
        else
            varargin(end) = [];
        end
    else
        % no parameter struct. so initialize an empty one.
        optParamStruct = struct();
    end
    nVarargin = numel(varargin); %added this in case varagin{1} is the parameter structure
    % now, loop through all parameter/value pairs.
    for i = 1:2:nVarargin
        cparam = varargin{i};
        if ~ischar(cparam)
            error('Parameters have to be supplied as ''parameterName''/Value pairs');
        end
        % the param struct overrides the only use the parameter if it is
        % not a field of the parameter struct.
        if ~isfield(optParamStruct,cparam)
            try
                optParamStruct.(cparam) = varargin{i+1};
            catch
                error('All parameters have to be valid matlab field names. %s is not a valid field name',cparam);
            end
        else
            warning('Duplicate parameter %s, both supplied as a field name and a parameter/value pair!',cparam);
        end
    end
else
    % no optional parameters.
    optParamStruct = struct();
end

% set up the cobra parameters
problemTypeParams = struct();

for i = 1:numel(defaultParams(:,1))
    % if the field is part of the optional parameters (i.e. explicitly provided) use it.
    if isfield(optParamStruct,defaultParams{i,1})
        problemTypeParams.(defaultParams{i,1}) = optParamStruct.(defaultParams{i,1});
        % and remove the field from the struct for the solver specific parameters.
        optParamStruct = rmfield(optParamStruct,defaultParams{i,1});
    else
        % otherwise use the default parameter
        problemTypeParams.(defaultParams{i,1}) = defaultParams{i,2};
    end
end

% assign all remaining parameters to the solver parameter struct.
solverParams = optParamStruct;

author: opencobra

repo: cobratoolbox


text:
function s = sudoku_line_row_swap ( box_row, line_row1, line_row2, s )

%*****************************************************************************80
%
%% SUDOKU_LINE_ROW_SWAP swaps two horizontal lines that share the same box row.
%
%  Discussion:
%
%    Two rows of a Sudoku may be shifted without changing the basic
%    structure, as long as the rows are in the same box row.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer BOX_ROW, the index of the box row, between 1 and 3.
%
%    Input, integer LINE_ROW1, LINE_ROW2, the indices of the two rows
%    which are to be interchanged.  These should be between 1 and 3.
%
%    Input, integer S(9,9), the Sudoku to be shuffled.
%
%    Output, integer S(9,9), the shuffled Sudoku.
%
  if ( 1 <= box_row & box_row <= 3 )
    if ( 1 <= line_row1 & line_row1 <= 3 )
      if ( 1 <= line_row2 & line_row2 <= 3 )

        i1 = 3 * ( box_row - 1 ) + line_row1;
        i2 = 3 * ( box_row - 1 ) + line_row2;

        t(1, 1:9) = s(i1,1:9);
        s(i1,1:9) = s(i2,1:9);
        s(i2,1:9) = t(1, 1:9);

      end
    end
  end

  return
end