PK J40 test_esn_8.mfunction [output, stateCollection] = test_esn_8(input, voteHistory, esn, nForgetInput) % test_esn_8 runs a trained ESN on a particular input % output has size [size(input,1) esn.d esn.p] % output contains the values of the output units *before* sigmoiding, with % reversing teacherShift and teacherScaling stateCollection = compute_statematrix_noTeacher_8(input, voteHistory, esn, nForgetInput) ; for n = 1:esn.d output(:,:,n) = feval(esn.dataActivationFunction, ... stateCollection * squeeze(esn.W_out(:,:,n))'); output(:,:,n) = ... output(:,:,n) - repmat(esn.teacherShift(:,n)', size(stateCollection,1),1); output(:,:,n) = ... output(:,:,n)./ repmat(esn.teacherScaling(:,n)', size(stateCollection,1),1); end PK q53e` ` compute_error_8.mfunction nrmse = compute_error_8(estimatedOutput, correctOutput, esn) % remark: estimatedOutput has size % (lengthTrainigData_minus_washout, dataSize, nOuts) % nrmse has size (dataSize, nOuts) % if length(correctOutput) > length(estimateOutput), the first % elements from correctOutput are deleted nEstimatedPoints = size(estimatedOutput, 1) ; nForgetPoints = size(correctOutput,1) - nEstimatedPoints ; correctOutputShort = correctOutput(nForgetPoints+1:end,:,:) ; mse = squeeze(mean(((estimatedOutput - correctOutputShort).^2), 1)) ; % size (dataSize, nOuts) variance = zeros(size(estimatedOutput,2), size(estimatedOutput,3)); % size (dataSize, nOuts) for d = 1:size(estimatedOutput,3) for p = 1:size(estimatedOutput,2) variance(p,d) = var(correctOutputShort(:,p,d)); end end nrmse = (sqrt(mse./variance)) ; PK v5< compute_statematrix_8.mfunction stateCollectMat = compute_statematrix_8(sampleInput, sampleOutput, voteHistory, ... esn, nForgetInput, trainNoiseLevel) nSampleInput = size(sampleInput,1); internalState = zeros(esn.reservoirSize, 1); totalstate = zeros(esn.reservoirSize + esn.p, 1); netOut = zeros(esn.p, esn.d); stateCollectMat = zeros(nSampleInput-nForgetInput, esn.reservoirSize + esn.p) ; W = esn.W * esn.spectralRadius; for iInput = 1:nSampleInput externalIn = sampleInput(iInput,:)'; % is column vector in = esn.inputScaling .* externalIn + esn.inputShift; voteState = [in netOut] * voteHistory(iInput,:)'; voteState = esn.inputWeightScaling .* voteState + esn.bias; totalstate(esn.reservoirSize+1:end,1) = voteState; internalState = feval(esn.reservoirActivationFunction, ... [W, esn.W_in] * totalstate ); internalState = internalState + trainNoiseLevel * (rand(esn.reservoirSize,1) - 0.5); % force teacher output teacher = squeeze(sampleOutput(iInput, :,:)); teacher = esn.teacherScaling .* teacher + esn.teacherShift; netOut = feval(esn.dataActivationFunction, teacher); totalstate(1:esn.reservoirSize) = internalState; %collect state if (iInput > nForgetInput) collectIndex = iInput - nForgetInput; stateCollectMat(collectIndex,:) = totalstate'; %fill a row end endPK 5P3 ! compute_statematrix_noTeacher_8.mfunction stateCollectMat = compute_statematrix_noTeacher_8(sampleInput, voteHistory, ... esn, nForgetInput) nSampleInput = size(sampleInput,1); internalState = zeros(esn.reservoirSize, 1); totalstate = zeros(esn.reservoirSize + esn.p, 1); netOut = zeros(esn.p, esn.d); stateCollectMat = zeros(nSampleInput-nForgetInput, esn.reservoirSize + esn.p) ; intWeights = esn.W * esn.spectralRadius; for iInput = 1:nSampleInput externalIn = sampleInput(iInput,:)'; % is column vector in = esn.inputScaling .* externalIn + esn.inputShift; voteState = [in netOut] * voteHistory(iInput,:)'; voteState = esn.inputWeightScaling .* voteState + esn.bias; totalstate(esn.reservoirSize+1:end,1) = voteState; internalState = feval(esn.reservoirActivationFunction, ... [intWeights, esn.W_in] * totalstate ); totalstate(1:esn.reservoirSize,1) = internalState; netOut = zeros(esn.p, esn.d); for n = 1:esn.d netOut(:,n) = (esn.W_out(:,:,n) * totalstate); end netOut = feval(esn.dataActivationFunction, netOut); %collect state if (iInput > nForgetInput) collectIndex = iInput - nForgetInput; stateCollectMat(collectIndex,:) = totalstate'; %fill a row end endPK WJ4 compute_teacher_8.mfunction teachCollectMat = compute_teacher_8(sampleOutput, esn, nForgetOutput) nSampleOutput = size(sampleOutput,1) - nForgetOutput ; teachCollectMat = zeros(nSampleOutput, esn.p, esn.d) ; %delete the first nForgetOutput elements from sampleOutput teachCollectMat = sampleOutput(nForgetOutput+1:end,:,:) ; % scaling and shifting for i = 1:nSampleOutput teachCollectMat(i,:,:) = feval(esn.inverseDataActivationFunction,... squeeze(teachCollectMat(i,:,:)) .* esn.teacherScaling + esn.teacherShift); end PK ˊ5AE E demo_8.m% Demo script for the IUB Technical Report Nr 3, Herbert Jaeger: % "Generating exponentially many periodic attractors with linearly growing Echo State % Networks". % % Herbert Jaeger, August 2006 % % Notice: International patents pending for the Echo State Network architecture % and learning method, claimed by Fraunhofer IAIS % http://www.iais.fraunhofer.de/ for most of Europe, USA, Canada, % Australia, Japan. p = 5; % nr of pitch bins d = 10; % nr of delays % 1. training ESN % 1.1 generate ESN reservoirSize = 100; % teacherScaling must be same as inputScaling, teacherShift same as % inputShift (except array sizes) esn = generate_esn_8(d, reservoirSize, p,'spectralRadius',0.9,... 'inputWeightScaling', 1.0 * ones(p,1),'bias', 0.0 * ones(p,1),... 'teacherScaling', 0.9 * ones(p, d), 'teacherShift', 0.05 * ones(p, d),... 'inputScaling', 0.9 * ones(p,1), 'inputShift', 0.05 * ones(p,1), ... 'reservoirActivationFunction', 'identity_8', ... 'dataActivationFunction', 'sigmoid01', ... 'inverseDataActivationFunction', 'sigmoid01_inv', ... 'methodTrain', 'wiener_hopf_8') ; % 1.2 prepare training data sampleLength = 1000; rawSignal = rand(sampleLength + d,1); sampleInput = rawToUnitCoded(rawSignal, p); sampleInput = cleanPitches(sampleInput); % order p x d (pitch times delay) outputs in d many portions each of length % p. Shortest delay come first sampleOutput = zeros(sampleLength, p, d); for i = 1:d sampleOutput(:,:,d-i+1) = sampleInput(i:end - d + i-1,:); end sampleInput = sampleInput(d + 1: end, :); [sampleInputTrain, sampleInputTest] = split_train_test_8(sampleInput, 0.75); [sampleOutputTrain, sampleOutputTest] = split_train_test_8(sampleOutput, 0.75); % 1.3 train the ESN nDiscard = 200; trainVotes = [ones(size(sampleInputTrain,1),1) zeros(size(sampleInputTrain,1), d)]; trainNoiseLevel = 0.0001; [esn.W_out stateMatrix teachMatrix] = ... train_esn_8(sampleInputTrain, sampleOutputTrain, trainVotes, esn, nDiscard, trainNoiseLevel) ; % 1.4 diagnostics for training figure(1); titlestring = 'four internal states during learning'; plot_states(stateMatrix,[1 2 3 4], 200, titlestring) ; figure(2); plot(squeeze(mean(abs(esn.W_out),2))); title('mean abs output weights'); estimatedTrainOutput = zeros(size(stateMatrix,1), p, d); for n = 1:esn.d estimatedTrainOutput(:,:,n) = feval(esn.dataActivationFunction, ... stateMatrix * squeeze(esn.W_out(:,:,n))'); estimatedTrainOutput(:,:,n) = ... estimatedTrainOutput(:,:,n) - repmat(esn.teacherShift(:,n)', size(stateMatrix,1),1); estimatedTrainOutput(:,:,n) = ... estimatedTrainOutput(:,:,n)./ repmat(esn.teacherScaling(:,n)', size(stateMatrix,1),1); end figure(3); clf; outputs_plotted = [1 1; 2 3; 3 5; 4 7; 5 10]; % select some pairs from p x d nplots = size(outputs_plotted, 1); lengthPlot = 200; for i = 1:nplots subplot(nplots, 1, i); plot(estimatedTrainOutput(1:lengthPlot, outputs_plotted(i,1), outputs_plotted(i,2))); hold on; plot(sampleOutputTrain(nDiscard+1:nDiscard+lengthPlot+1, ... outputs_plotted(i,1), outputs_plotted(i,2)), 'g:'); hold off; if i == 1 title('trained ESN output(b) vs. teacher(g:)'); end end trainNRMSE = compute_error_8(estimatedTrainOutput, sampleOutputTrain, esn); % size (p, d) figure(4); clf; plot(trainNRMSE); title('training NRMSE'); %compute NMSE testing error voteHistory = [ones(size(sampleInputTest,1),1) zeros(size(sampleInputTest,1), d)]; [testOutput, stateCollection] = test_esn_8(sampleInputTest, voteHistory, esn, nDiscard) ; testNRMSE = compute_error_8(testOutput, sampleOutputTest, esn); figure(5); clf; plot(testNRMSE); title('testing NRMSE'); % use delay-line ESN to continue rhythmic trigger pattern initN = 20; periodN = 7; repeatN = 10; triggerInputRaw = rand(initN + periodN * repeatN, 1) ; periodPattern = rand(periodN ,1); for r = 1:repeatN triggerInputRaw(initN + (r-1)*periodN + 1 : initN + r * periodN, 1) = periodPattern; end triggerInput = rawToUnitCoded(triggerInputRaw, esn.p); triggerInput = cleanPitches(triggerInput); triggerInputRaw = unitCodedToRaw(triggerInput); zeroConfidenceThreshold = 0.3; votingGrowthConstant = 4; votingDecay = 0.2; errDecay = 0.4; errGrowthConstant = 4; startAfterRepeatsN = 2; noiseLevel = 0.01; netOutCompletePL = run_melody_8(esn, triggerInput, triggerInputRaw, ... initN, periodN, repeatN, startAfterRepeatsN, ... zeroConfidenceThreshold, votingGrowthConstant, votingDecay, errDecay, ... errGrowthConstant, noiseLevel); PK L4!2 2 generate_esn_8.mfunction esn = generate_esn_8(d, reservoirSize, p, varargin) % generate an untrained "raw" ESN with p input neurons, p * d output % neurons. Idea: p neurons code pitch, d is number of trained delays. clear esn; esn.reservoirSize = reservoirSize; esn.d = d; esn.p = p; % create internal weight matrix with spectral radius 1 and connectivity % such that on average each unit is connected to 10 others connectivity = min([10/reservoirSize 1]); esn.W = generate_internal_weights(reservoirSize, connectivity); % input weight matrix has weight vectors per input unit in colums esn.W_in = 2.0 * rand(reservoirSize, p)- 1.0; % initialize output weights esn.W_out = zeros(p, reservoirSize + p, d); %init default parameters esn.spectralRadius = 0.8; esn.inputWeightScaling = ones(p,1); esn.bias = zeros(p,1); esn.teacherScaling = ones(p, d); esn.teacherShift = zeros(p, d); esn.inputScaling = ones(p, 1); esn.inputShift = zeros(p, 1); esn.noiseLevel = 0.0 ; esn.reservoirActivationFunction = 'tanh'; esn.dataActivationFunction = 'sigmoid01' ; esn.inverseDataActivationFunction = 'sigmoid01_inv' ; esn.methodTrain = 'pseudoinverse_8' ; args = varargin; nargs= length(args); for i = 1:2:nargs switch args{i}, case 'spectralRadius', esn.spectralRadius = args{i+1} ; case 'inputWeightScaling', esn.inputWeightScaling = args{i+1} ; case 'bias', esn.bias = args{i+1} ; case 'teacherScaling', esn.teacherScaling = args{i+1} ; case 'teacherShift', esn.teacherShift = args{i+1} ; case 'inputScaling', esn.inputScaling = args{i+1} ; case 'inputShift', esn.inputShift = args{i+1} ; case 'noiseLevel', esn.noiseLevel = args{i+1} ; case 'reservoirActivationFunction', esn.reservoirActivationFunction = args{i+1}; case 'dataActivationFunction', esn.dataActivationFunction = args{i+1}; case 'inverseDataActivationFunction', esn.inverseDataActivationFunction = args{i+1}; case 'methodTrain', esn.methodTrain = args{i+1} ; otherwise error('the option does not exist'); end end PK H 49&H generate_internal_weights.mfunction internalWeights = generate_internal_weights(nInternalUnits, ... connectivity) % GENERATE_INTERNAL_WEIGHTS creates a random reservoir for an ESN % % inputs: % nInternalUnits = the number of internal units in the ESN % connectivity \in [0,1], says how many weights should be non-zero % % output: % internalWeights = matrix of size nInternalUnits x nInternalUnits % internalWeights(i,j) = value of weight(synapse) from unit i to unit j % internalWeights(i,j) might be different from internalWeights(j,i) % % Version 1.0, April 30, 2006 % (c) Dan Popovici % % supress any diagnostic messages when computing the largest eigenvalue opts.disp = 0 ; succes = 0 ; while succes == 0 % following block might fail, thus we repeat until we obtain a valid % internalWeights matrix try, internalWeights = sprand(nInternalUnits, nInternalUnits, connectivity); internalWeights = spfun(@minusPoint5,internalWeights); maxVal = max(abs(eigs(internalWeights,1,'LM',opts))); internalWeights = internalWeights/maxVal; succes = 1 ; catch, succes = 0 ; end end PK 5q" " identity_8.mfunction y = identity_8(x) y = x;PK #}/=_) ) minusPoint5.mfunction x = minusPoint5(y) x = y - 0.5;PK }/8ҁ] ] myeigs.mfunction varargout = eigs(varargin) %EIGS Find a few eigenvalues and eigenvectors of a matrix using ARPACK. % D = EIGS(A) returns a vector of A's 6 largest magnitude eigenvalues. % A must be square and should be large and sparse. % % [V,D] = EIGS(A) returns a diagonal matrix D of A's 6 largest magnitude % eigenvalues and a matrix V whose columns are the corresponding eigenvectors. % % [V,D,FLAG] = EIGS(A) also returns a convergence flag. If FLAG is 0 % then all the eigenvalues converged; otherwise not all converged. % % EIGS(A,B) solves the generalized eigenvalue problem A*V == B*V*D. B must % be symmetric (or Hermitian) positive definite and the same size as A. % EIGS(A,[],...) indicates the standard eigenvalue problem A*V == V*D. % % EIGS(A,K) and EIGS(A,B,K) return the K largest magnitude eigenvalues. % % EIGS(A,K,SIGMA) and EIGS(A,B,K,SIGMA) return K eigenvalues based on SIGMA: % 'LM' or 'SM' - Largest or Smallest Magnitude % For real symmetric problems, SIGMA may also be: % 'LA' or 'SA' - Largest or Smallest Algebraic % 'BE' - Both Ends, one more from high end if K is odd % For nonsymmetric and complex problems, SIGMA may also be: % 'LR' or 'SR' - Largest or Smallest Real part % 'LI' or 'SI' - Largest or Smallest Imaginary part % If SIGMA is a real or complex scalar including 0, EIGS finds the eigenvalues % closest to SIGMA. For scalar SIGMA, and also when SIGMA = 'SM' which uses % the same algorithm as SIGMA = 0, B need only be symmetric (or Hermitian) % positive semi-definite since it is not Cholesky factored as in the other cases. % % EIGS(A,K,SIGMA,OPTS) and EIGS(A,B,K,SIGMA,OPTS) specify options: % OPTS.issym: symmetry of A or A-SIGMA*B represented by AFUN [{0} | 1] % OPTS.isreal: complexity of A or A-SIGMA*B represented by AFUN [0 | {1}] % OPTS.tol: convergence: Ritz estimate residual <= tol*NORM(A) [scalar | {eps}] % OPTS.maxit: maximum number of iterations [integer | {300}] % OPTS.p: number of Lanczos vectors: K+1
3) error('Too many output arguments.') end % Process inputs and do error-checking if isa(varargin{1},'double') A = varargin{1}; Amatrix = 1; else A = fcnchk(varargin{1}); Amatrix = 0; end isrealprob = 1; % isrealprob = isreal(A) & isreal(B) & isreal(sigma) if Amatrix isrealprob = isreal(A); end issymA = 0; if Amatrix issymA = isequal(A,A'); end if Amatrix [m,n] = size(A); if (m ~= n) error('A must be a square matrix or a function which computes A*x.') end else n = varargin{2}; nstr = 'Size of problem, ''n'', must be a positive integer.'; if ~isequal(size(n),[1,1]) | ~isreal(n) error(nstr) end if (round(n) ~= n) warning('MATLAB:eigs:NonIntegerSize',['%s\n ' ... 'Rounding input size.'],nstr) n = round(n); end if issparse(n) n = full(n); end end Bnotthere = 0; Bstr = sprintf(['Generalized matrix B must be the same size as A and' ... ' either a symmetric positive (semi-)definite matrix or' ... ' its Cholesky factor.']); if (nargin < (3-Amatrix-Bnotthere)) B = []; Bnotthere = 1; else Bk = varargin{3-Amatrix-Bnotthere}; if isempty(Bk) % allow eigs(A,[],k,sigma,opts); B = Bk; else if isequal(size(Bk),[1,1]) & (n ~= 1) B = []; k = Bk; Bnotthere = 1; else % eigs(9,8,...) assumes A=9, B=8, ... NOT A=9, k=8, ... B = Bk; if ~isa(B,'double') | ~isequal(size(B),[n,n]) error(Bstr) end isrealprob = isrealprob & isreal(B); end end end if Amatrix & ((nargin - ~Bnotthere)>4) error('Too many inputs.') end if (nargin < (4-Amatrix-Bnotthere)) k = min(n,6); else k = varargin{4-Amatrix-Bnotthere}; end kstr = ['Number of eigenvalues requested, k, must be a' ... ' positive integer <= n.']; if ~isa(k,'double') | ~isequal(size(k),[1,1]) | ~isreal(k) | (k>n) error(kstr) end if issparse(k) k = full(k); end if (round(k) ~= k) warning('MATLAB:eigs:NonIntegerEigQty',['%s\n ' ... 'Rounding number of eigenvalues.'],kstr) k = round(k); end whchstr = sprintf(['Eigenvalue range sigma must be a valid 2-element string.']); if (nargin < (5-Amatrix-Bnotthere)) % default: eigs('LM') => ARPACK which='LM', sigma=0 eigs_sigma = 'LM'; whch = 'LM'; sigma = 0; else eigs_sigma = varargin{5-Amatrix-Bnotthere}; if isstr(eigs_sigma) % eigs(string) => ARPACK which=string, sigma=0 if ~isequal(size(eigs_sigma),[1,2]) whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ... 'LM','SM','LA','SA','BE')]; whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ... 'LM','SM','LR','SR','LI','SI')]; error(whchstr) end eigs_sigma = upper(eigs_sigma); if isequal(eigs_sigma,'SM') % eigs('SM') => ARPACK which='LM', sigma=0 whch = 'LM'; else % eigs(string), where string~='SM' => ARPACK which=string, sigma=0 whch = eigs_sigma; end sigma = 0; else % eigs(scalar) => ARPACK which='LM', sigma=scalar if ~isa(eigs_sigma,'double') | ~isequal(size(eigs_sigma),[1,1]) error('Eigenvalue shift sigma must be a scalar.') end sigma = eigs_sigma; if issparse(sigma) sigma = full(sigma); end isrealprob = isrealprob & isreal(sigma); whch = 'LM'; end end tol = eps; % ARPACK's minimum tolerance is eps/2 (DLAMCH's EPS) maxit = []; p = []; info = int32(0); % use a random starting vector display = 0; cholB = 0; if (nargin >= (6-Amatrix-Bnotthere)) opts = varargin{6-Amatrix-Bnotthere}; if ~isa(opts,'struct') error('Options argument must be a structure.') end if isfield(opts,'issym') & ~Amatrix issymA = opts.issym; if (issymA ~= 0) & (issymA ~= 1) error('opts.issym must be 0 or 1.') end end if isfield(opts,'isreal') & ~Amatrix if (opts.isreal ~= 0) & (opts.isreal ~= 1) error('opts.isreal must be 0 or 1.') end isrealprob = isrealprob & opts.isreal; end if ~isempty(B) & (isfield(opts,'cholB') | isfield(opts,'permB')) if isfield(opts,'cholB') cholB = opts.cholB; if (cholB ~= 0) & (cholB ~= 1) error('opts.cholB must be 0 or 1.') end if isfield(opts,'permB') if issparse(B) & cholB permB = opts.permB; if ~isequal(sort(permB),(1:n)) & ... ~isequal(sort(permB),(1:n)') error('opts.permB must be a permutation of 1:n.') end else warning('MATLAB:eigs:IgnoredOptionPermB', ... ['Ignoring opts.permB since B is not its sparse' ... ' Cholesky factor.']) end else permB = 1:n; end end end if isfield(opts,'tol') if ~isequal(size(opts.tol),[1,1]) | ~isreal(opts.tol) | (opts.tol<=0) error(['Convergence tolerance opts.tol must be a strictly' ... ' positive real scalar.']) else tol = full(opts.tol); end end if isfield(opts,'p') p = opts.p; pstr = ['Number of basis vectors opts.p must be a positive' ... ' integer <= n.']; if ~isequal(size(p),[1,1]) | ~isreal(p) | (p<=0) | (p>n) error(pstr) end if issparse(p) p = full(p); end if (round(p) ~= p) warning('MATLAB:eigs:NonIntegerVecQty',['%s\n ' ... 'Rounding number of basis vectors.'],pstr) p = round(p); end end if isfield(opts,'maxit') maxit = opts.maxit; str = ['Maximum number of iterations opts.maxit must be' ... ' a positive integer.']; if ~isequal(size(maxit),[1,1]) | ~isreal(maxit) | (maxit<=0) error(str) end if issparse(maxit) maxit = full(maxit); end if (round(maxit) ~= maxit) warning('MATLAB:eigs:NonIntegerIterationQty',['%s\n ' ... 'Rounding number of iterations.'],str) maxit = round(maxit); end end if isfield(opts,'v0') if ~isequal(size(opts.v0),[n,1]) error('Start vector opts.v0 must be n-by-1.') end if isrealprob if ~isreal(opts.v0) error(['Start vector opts.v0 must be real for real problems.']) end resid = full(opts.v0); else resid(1:2:(2*n-1),1) = full(real(opts.v0)); resid(2:2:2*n,1) = full(imag(opts.v0)); end info = int32(1); % use resid as starting vector end if isfield(opts,'disp') display = opts.disp; dispstr = 'Diagnostic level opts.disp must be an integer.'; if (~isequal(size(display),[1,1])) | (~isreal(display)) | (display<0) error(dispstr) end if (round(display) ~= display) warning('MATLAB:eigs:NonIntegerDiagnosticLevel', ... '%s\n Rounding diagnostic level.',dispstr) display = round(display); end end if isfield(opts,'cheb') warning('MATLAB:eigs:ObsoleteOptionCheb', ... ['Ignoring polynomial acceleration opts.cheb' ... ' (no longer an option).']); end if isfield(opts,'stagtol') warning('MATLAB:eigs:ObsoleteOptionStagtol', ... ['Ignoring stagnation tolerance opts.stagtol' ... ' (no longer an option).']); end end % Now we know issymA, isrealprob, cholB, and permB if isempty(p) if isrealprob & ~issymA p = min(max(2*k+1,20),n); else p = min(max(2*k,20),n); end end if isempty(maxit) maxit = max(300,ceil(2*n/max(p,1))); end if (info == int32(0)) if isrealprob resid = zeros(n,1); else resid = zeros(2*n,1); end end if ~isempty(B) % B must be symmetric (Hermitian) positive (semi-)definite if cholB if ~isequal(triu(B),B) error(Bstr) end else if ~isequal(B,B') error(Bstr) end end end useeig = 0; if isrealprob & issymA knstr = sprintf(['For real symmetric problems, must have number' ... ' of eigenvalues k < n.\n']); else knstr = sprintf(['For nonsymmetric and complex problems, must have' ... ' number of eigenvalues k < n-1.\n']); end if isempty(B) knstr = [knstr sprintf(['Using eig(full(A)) instead.'])]; else knstr = [knstr sprintf(['Using eig(full(A),full(B)) instead.'])]; end if (k == 0) useeig = 1; end if isrealprob & issymA if (k > n-1) if (n >= 6) warning('MATLAB:eigs:TooManyRequestedEigsForRealSym', ... '%s',knstr) end useeig = 1; end else if (k > n-2) if (n >= 7) warning('MATLAB:eigs:TooManyRequestedEigsForComplexNonsym', ... '%s',knstr) end useeig = 1; end end if isrealprob & issymA if ~isreal(sigma) error(['For real symmetric problems, eigenvalue shift sigma must' ... ' be real.']) end else if ~isrealprob & issymA & ~isreal(sigma) warning('MATLAB:eigs:ComplexShiftForRealProblem', ... ['Complex eigenvalue shift sigma on a Hermitian problem' ... ' (all real eigenvalues).']) end end if isrealprob & issymA if strcmp(whch,'LR') whch = 'LA'; warning('MATLAB:eigs:SigmaChangedToLA', ... ['For real symmetric problems, sigma value ''LR''' ... ' (Largest Real) is now ''LA'' (Largest Algebraic).']) end if strcmp(whch,'SR') whch = 'SA'; warning('MATLAB:eigs:SigmaChangedToSA', ... ['For real symmetric problems, sigma value ''SR''' ... ' (Smallest Real) is now ''SA'' (Smallest Algebraic).']) end if ~ismember(whch,{'LM', 'SM', 'LA', 'SA', 'BE'}) whchstr = [whchstr sprintf(['\nFor real symmetric A, the choices are ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.'], ... 'LM','SM','LA','SA','BE')]; error(whchstr) end else if strcmp(whch,'BE') warning('MATLAB:eigs:SigmaChangedToLM', ... ['Sigma value ''BE'' is now only available for real' ... ' symmetric problems. Computing ''LM'' eigenvalues instead.']) whch = 'LM'; end if ~ismember(whch,{'LM', 'SM', 'LR', 'SR', 'LI', 'SI'}) whchstr = [whchstr sprintf(['\nFor non-symmetric or complex A, the choices are ''%s'', ''%s'', ''%s'', ''%s'', ''%s'' or ''%s''.\n'], ... 'LM','SM','LR','SR','LI','SI')]; error(whchstr) end end % Now have enough information to do early return on cases eigs does not handle if useeig if (nargout <= 1) varargout{1} = eigs2(A,n,B,k,whch,sigma,cholB, ... varargin{7-Amatrix-Bnotthere:end}); else [varargout{1},varargout{2}] = eigs2(A,n,B,k,whch,sigma,cholB, ... varargin{7-Amatrix-Bnotthere:end}); end if (nargout == 3) varargout{3} = 0; % flag indicates "convergence" end return end if isrealprob & ~issymA sigmar = real(sigma); sigmai = imag(sigma); end if isrealprob & issymA if (p <= k) error(['For real symmetric problems, must have number of' ... ' basis vectors opts.p > k.']) end else if (p <= k+1) error(['For nonsymmetric and complex problems, must have number of' ... ' basis vectors opts.p > k+1.']) end end if isequal(whch,'LM') & ~isequal(eigs_sigma,'LM') % A*x = lambda*M*x, M symmetric (positive) semi-definite % => OP = inv(A - sigma*M)*M and B = M % => shift-and-invert mode mode = 3; elseif isempty(B) % A*x = lambda*x % => OP = A and B = I mode = 1; else % B is not empty % Do not use mode=2. % Use mode = 1 with OP = R'\(A*(R\x)) and B = I % where R is B's upper triangular Cholesky factor: B = R'*R. % Finally, V = R\V returns the actual generalized eigenvectors of A and B. mode = 1; end if cholB pB = 0; RB = B; RBT = B'; end if (mode == 3) & Amatrix % need lu(A-sigma*B) if issparse(A) & (isempty(B) | issparse(B)) if isempty(B) AsB = A - sigma * speye(n); else if cholB AsB = A - sigma * RBT * RB; else AsB = A - sigma * B; end end [L,U,P,Q] = lu(AsB); [perm,dummy] = find(Q); else if isempty(B) AsB = A - sigma * eye(n); else if cholB AsB = A - sigma * RBT * RB; else AsB = A - sigma * B; end end [L,U,P] = lu(AsB); end dU = diag(U); rcondestU = full(min(abs(dU)) / max(abs(dU))); if (rcondestU < eps) if isempty(B) ds = sprintf(['(A-sigma*I) has small reciprocal condition' ... ' estimate: %f\n'],rcondestU); else ds = sprintf(['(A-sigma*B) has small reciprocal condition' ... ' estimate: %f\n'],rcondestU); end ds = [ds sprintf(['indicating that sigma is near an exact' ... ' eigenvalue. The\nalgorithm may not converge unless' ... ' you try a new value for sigma.\n'])]; warning('MATLAB:eigs:SigmaNearExactEig',ds) end end if (mode == 1) & ~isempty(B) & ~cholB % need chol(B) if issparse(B) permB = symamd(B); [RB,pB] = chol(B(permB,permB)); else [RB,pB] = chol(B); end if (pB == 0) RBT = RB'; else error(Bstr) end end % Allocate outputs and ARPACK work variables if isrealprob if issymA % real and symmetric prefix = 'ds'; v = zeros(n,p); ldv = int32(size(v,1)); ipntr = int32(zeros(15,1)); workd = zeros(n,3); lworkl = p*(p+8); workl = zeros(lworkl,1); lworkl = int32(lworkl); d = zeros(k,1); else % real but not symmetric prefix = 'dn'; v = zeros(n,p); ldv = int32(size(v,1)); ipntr = int32(zeros(15,1)); workd = zeros(n,3); lworkl = 3*p*(p+2); workl = zeros(lworkl,1); lworkl = int32(lworkl); workev = zeros(3*p,1); d = zeros(k+1,1); di = zeros(k+1,1); end else % complex prefix = 'zn'; zv = zeros(2*n*p,1); ldv = int32(n); ipntr = int32(zeros(15,1)); workd = complex(zeros(n,3)); zworkd = zeros(2*prod(size(workd)),1); lworkl = 3*p^2+5*p; workl = zeros(2*lworkl,1); lworkl = int32(lworkl); workev = zeros(2*2*p,1); zd = zeros(2*(k+1),1); rwork = zeros(p,1); end ido = int32(0); % reverse communication parameter if isempty(B) | (~isempty(B) & (mode == 1)) bmat = 'I'; % standard eigenvalue problem else bmat = 'G'; % generalized eigenvalue problem end nev = int32(k); % number of eigenvalues requested ncv = int32(p); % number of Lanczos vectors iparam = int32(zeros(11,1)); iparam([1 3 7]) = int32([1 maxit mode]); select = int32(zeros(p,1)); cputms(1) = cputime - t0; % end timing pre-processing iter = 0; ariter = 0; while (ido ~= 99) t0 = cputime; % start timing ARPACK calls **aupd if isrealprob arpackc( [prefix 'aupd'], ido, ... bmat, int32(n), whch, nev, tol, resid, ncv, ... v, ldv, iparam, ipntr, workd, workl, lworkl, info); else zworkd(1:2:end-1) = real(workd); zworkd(2:2:end) = imag(workd); arpackc( 'znaupd', ido, ... bmat, int32(n), whch, nev, tol, resid, ncv, zv, ... ldv, iparam, ipntr, zworkd, workl, lworkl, ... rwork, info ); workd = reshape(complex(zworkd(1:2:end-1),zworkd(2:2:end)),[n,3]); end if (info < 0) es = sprintf('Error with ARPACK routine %saupd: info = %d',... prefix,full(info)); error(es) end cputms(2) = cputms(2) + (cputime-t0); % end timing ARPACK calls **aupd t0 = cputime; % start timing MATLAB OP(X) % Compute which columns of workd ipntr references [row,col1] = ind2sub([n,3],double(ipntr(1))); if (row ~= 1) str = sprintf(['ipntr(1)=%d does not refer to the start of a' ... ' column of the %d-by-3 array workd.'],full(ipntr(1)),n); error(str) end [row,col2] = ind2sub([n,3],double(ipntr(2))); if (row ~= 1) str = sprintf(['ipntr(2)=%d does not refer to the start of a' ... ' column of the %d-by-3 array workd.'],full(ipntr(2)),n); error(str) end if ~isempty(B) & (mode == 3) & (ido == 1) [row,col3] = ind2sub([n,3],double(ipntr(3))); if (row ~= 1) str = sprintf(['ipntr(3)=%d does not refer to the start of a' ... ' column of the %d-by-3 array workd.'],full(ipntr(3)),n); error(str) end end switch (ido) case {-1,1} if Amatrix if (mode == 1) if isempty(B) % OP = A*x workd(:,col2) = A * workd(:,col1); else % OP = R'\(A*(R\x)) if issparse(B) workd(permB,col2) = RB \ workd(:,col1); workd(:,col2) = A * workd(:,col2); workd(:,col2) = RBT \ workd(permB,col2); else workd(:,col2) = RBT \ (A * (RB \ workd(:,col1))); end end elseif (mode == 3) if isempty(B) if issparse(A) workd(perm,col2) = U \ (L \ (P * workd(:,col1))); else workd(:,col2) = U \ (L \ (P * workd(:,col1))); end else % B is not empty if (ido == -1) if cholB workd(:,col2) = RBT * (RB * workd(:,col1)); else workd(:,col2) = B * workd(:,col1); end if issparse(A) & issparse(B) workd(perm,col2) = U \ (L \ (P * workd(:,col1))); else workd(:,col2) = U \ (L \ (P * workd(:,col1))); end elseif (ido == 1) if issparse(A) & issparse(B) workd(perm,col2) = U \ (L \ (P * workd(:,col3))); else workd(:,col2) = U \ (L \ (P * workd(:,col3))); end end end else % mode is not 1 or 3 error(['Unknown mode returned from ' prefix 'aupd.']) end else % A is not a matrix if (mode == 1) if isempty(B) % OP = A*x workd(:,col2) = feval(A,workd(:,col1), ... varargin{7-Amatrix-Bnotthere:end}); else % OP = R'\(A*(R\x)) if issparse(B) workd(permB,col2) = RB \ workd(:,col1); workd(:,col2) = feval(A,workd(:,col2),varargin{7-Amatrix-Bnotthere:end}); workd(:,col2) = RBT \ workd(permB,col2); else workd(:,col2) = RBT \ ... feval(A,(RB\workd(:,col1)), ... varargin{7-Amatrix-Bnotthere:end}); end end elseif (mode == 3) if isempty(B) workd(:,col2) = feval(A,workd(:,col1), ... varargin{7-Amatrix-Bnotthere:end}); else if (ido == -1) if cholB workd(:,col2) = RBT * (RB * workd(:,col1)); else workd(:,col2) = B * workd(:,col1); end workd(:,col2) = feval(A,workd(:,col2), ... varargin{7-Amatrix-Bnotthere:end}); elseif (ido == 1) workd(:,col2) = feval(A,workd(:,col3), ... varargin{7-Amatrix-Bnotthere:end}); end end else % mode is not 1 or 3 error(['Unknown mode returned from ' prefix 'aupd.']) end end % if Amatrix case 2 if (mode == 3) if cholB workd(:,col2) = RBT * (RB * workd(:,col1)); else workd(:,col2) = B * workd(:,col1); end else error(['Unknown mode returned from ' prefix 'aupd.']) end case 3 % setting iparam(1) = ishift = 1 ensures this never happens warning('MATLAB:eigs:WorklShiftsUnsupported', ... ['EIGS does not yet support computing the shifts in workl.' ... ' Setting reverse communication parameter to 99 and returning.']) ido = int32(99); case 99 otherwise error(['Unknown value of reverse communication parameter' ... ' returned from ' prefix 'aupd.']) end cputms(3) = cputms(3) + (cputime-t0); % end timing MATLAB OP(X) if display iter = double(ipntr(15)); if (iter > ariter) & (ido ~= 99) ariter = iter; % ds = sprintf(['Iteration %d: a few Ritz values of the' ... % ' %d-by-%d matrix:'],iter,p,p); % disp(ds) if isrealprob if issymA dispvec = [workl(double(ipntr(6))+(0:p-1))]; if isequal(whch,'BE') % roughly k Large eigenvalues and k Small eigenvalues disp(dispvec(max(end-2*k+1,1):end)) else % k eigenvalues disp(dispvec(max(end-k+1,1):end)) end else dispvec = [complex(workl(double(ipntr(6))+(0:p-1)), ... workl(double(ipntr(7))+(0:p-1)))]; % k+1 eigenvalues (keep complex conjugate pairs together) disp(dispvec(max(end-k,1):end)) end else dispvec = [complex(workl(2*double(ipntr(6))-1+(0:2:2*(p-1))), ... workl(2*double(ipntr(6))+(0:2:2*(p-1))))]; disp(dispvec(max(end-k+1,1):end)) end end end end % while (ido ~= 99) t0 = cputime; % start timing post-processing flag = 0; if (info < 0) es = sprintf('Error with ARPACK routine %saupd: info = %d',prefix,full(info)); error(es) else if (nargout >= 2) rvec = int32(1); % compute eigenvectors else rvec = int32(0); % do not compute eigenvectors end if isrealprob if issymA arpackc( 'dseupd', rvec, 'A', select, ... d, v, ldv, sigma, ... bmat, int32(n), whch, nev, tol, resid, ncv, ... v, ldv, iparam, ipntr, workd, workl, lworkl, info ); if isequal(whch,'LM') | isequal(whch,'LA') d = flipud(d); if (rvec == 1) v(:,1:k) = v(:,k:-1:1); end end if ((isequal(whch,'SM') | isequal(whch,'SA')) & (rvec == 0)) d = flipud(d); end else arpackc( 'dneupd', rvec, 'A', select, ... d, di, v, ldv, sigmar, sigmai, workev, ... bmat, int32(n), whch, nev, tol, resid, ncv, ... v, ldv, iparam, ipntr, workd, workl, lworkl, info ); d = complex(d,di); if rvec d(k+1) = []; else zind = find(d == 0); if isempty(zind) d = d(k+1:-1:2); else d(max(zind)) = []; d = flipud(d); end end end else zsigma = [real(sigma); imag(sigma)]; arpackc( 'zneupd', rvec, 'A', select, ... zd, zv, ldv, zsigma, workev, ... bmat, int32(n), whch, nev, tol, resid, ncv, zv, ... ldv, iparam, ipntr, zworkd, workl, lworkl, ... rwork, info ); if issymA d = zd(1:2:end-1); else d = complex(zd(1:2:end-1),zd(2:2:end)); end v = reshape(complex(zv(1:2:end-1),zv(2:2:end)),[n p]); end if (info ~= 0) es = ['Error with ARPACK routine ' prefix 'eupd: ']; switch double(info) case 2 ss = sum(select); if (ss < k) es = [es ... ' The logical variable select was only set with ' int2str(ss) ... ' 1''s instead of nconv=' int2str(double(iparam(5))) ... ' (k=' int2str(k) ').' ... ' Please contact the ARPACK authors at arpack@caam.rice.edu.']; else es = [es ... 'The LAPACK reordering routine ' prefix(1) ... 'trsen did not return all ' int2str(k) ' eigenvalues.'] end case 1 es = [es ... 'The Schur form could not be reordered by the LAPACK routine ' ... prefix(1) 'trsen.' ... ' Please contact the ARPACK authors at arpack@caam.rice.edu.']; case -14 es = [es prefix ... 'aupd did not find any eigenvalues to sufficient accuracy.']; otherwise es = [es sprintf('info = %d',full(info))]; end error(es) else nconv = double(iparam(5)); if (nconv == 0) if (nargout < 3) warning('MATLAB:eigs:NoEigsConverged', ... 'None of the %d requested eigenvalues converged.',k) else flag = 1; end elseif (nconv < k) if (nargout < 3) warning('MATLAB:eigs:NotAllEigsConverged', ... 'Only %d of the %d requested eigenvalues converged.',nconv,k) else flag = 1; end end end % if (info ~= 0) end % if (info < 0) if (issymA) | (~isrealprob) if (nargout <= 1) if isrealprob varargout{1} = d; else varargout{1} = d(k:-1:1,1); end else varargout{1} = v(:,1:k); varargout{2} = diag(d(1:k,1)); if (nargout >= 3) varargout{3} = flag; end end else if (nargout <= 1) varargout{1} = d; else cplxd = find(di ~= 0); % complex conjugate pairs of eigenvalues occur together cplxd = cplxd(1:2:end); v(:,[cplxd cplxd+1]) = [complex(v(:,cplxd),v(:,cplxd+1)) ... complex(v(:,cplxd),-v(:,cplxd+1))]; varargout{1} = v(:,1:k); varargout{2} = diag(d); if (nargout >= 3) varargout{3} = flag; end end end if (nargout >= 2) & (mode == 1) & ~isempty(B) if issparse(B) varargout{1}(permB,:) = RB \ varargout{1}; else varargout{1} = RB \ varargout{1}; end end cputms(4) = cputime-t0; % end timing post-processing cputms(5) = sum(cputms(1:4)); % total time if (display == 2) if (mode == 1) innerstr = sprintf(['Compute A*X:' ... ' %f\n'],cputms(3)); elseif (mode == 2) innerstr = sprintf(['Compute A*X and solve B*X=Y for X:' ... ' %f\n'],cputms(3)); elseif (mode == 3) if isempty(B) innerstr = sprintf(['Solve (A-SIGMA*I)*X=Y for X:' ... ' %f\n'],cputms(3)); else innerstr = sprintf(['Solve (A-SIGMA*B)*X=B*Y for X:' ... ' %f\n'],cputms(3)); end end if ((mode == 3) & (Amatrix)) if isempty(B) prepstr = sprintf(['Pre-processing, including lu(A-sigma*I):' ... ' %f\n'],cputms(1)); else prepstr = sprintf(['Pre-processing, including lu(A-sigma*B):' ... ' %f\n'],cputms(1)); end elseif ((mode == 2) & (~cholB)) prepstr = sprintf(['Pre-processing, including chol(B):' ... ' %f\n'],cputms(1)); else prepstr = sprintf(['Pre-processing:' ... ' %f\n'],cputms(1)); end sstr = sprintf(['***********CPU Timing Results in seconds***********']); ds = sprintf(['\n' sstr '\n' ... prepstr ... 'ARPACK''s %saupd: %f\n' ... innerstr ... 'Post-processing with ARPACK''s %seupd: %f\n' ... '***************************************************\n' ... 'Total: %f\n' ... sstr '\n'], ... prefix,cputms(2),prefix,cputms(4),cputms(5)); disp(ds) end PK o5Iع3 3 plot_states.mfunction plot_states(stateMatrix, states, nPoints, titlestring) %PLOT_STATES plots the internal states of the esn %example : plot_states(stateMatrix,[1 2 3 4],200) plots the first 200 %points from the traces of the first 4 units if (nargin < 3) nPoints = length(stateMatrix(:,1)) ; end nStates = length(states) ; xMax = ceil(sqrt(nStates)) ; yMax = ceil(nStates /xMax); for iPlot = 1 : nStates subplot(xMax,xMax,iPlot) ; plot(stateMatrix(1:nPoints,iPlot)); if iPlot == 1 title(titlestring); end end PK aL4K8 pseudoinverse_8.mfunction outputWeights = pseudoinverse_8(stateCollectMat, teachCollectMat, esn) pinvSCM = pinv(stateCollectMat); for i = 1:esn.d outputWeights(:,:,i) = (pinvSCM * teachCollectMat(:,:,i))' ; end PK 3e rawToUnitCoded.mfunction unitCodedSignal = rawToUnitCoded(rawSignal, nPitches) % transforms a 1-d stretchedSignal rawSignal with range [0,1] into % nPitches-dimensional signal % (each dim ranging in [0,1]) through nPitches many triangular, equally % spaced membership functions % % unitCodedSignal has size (length(rawSignal) , nPitches) l = length(rawSignal); unitCodedSignal = zeros(l, nPitches); stretchedSignal = rawSignal * (nPitches - 1); for i = 1:l for p = 1:nPitches unitCodedSignal(i, p) = unitTriangle(stretchedSignal(i,1) - p + 1); end end function membershipValue = unitTriangle(x) % computes the membership function of a unit triangle centered at zero if abs(x) > 1 membershipValue = 0; else membershipValue = 1 - abs(x); endPK 85Ko