ACHRSamplerParallelGeneral

PURPOSE ^

ACHRSamplerParallelGeneral Artificial Centering Hit-and-Run sampler with in place (memory) point

SYNOPSIS ^

function [sampleStruct] = ACHRSamplerParallelGeneral(sampleStruct,nLoops,stepsPerPoint, maxtime, proc, fdirectory)

DESCRIPTION ^

 ACHRSamplerParallelGeneral Artificial Centering Hit-and-Run sampler with in place (memory) point
 managmenet

 sampleStruct = ACHRSamplerParallelGeneral(sampleStruct,nLoops,stepsPerPoint)

INPUTS
 sampleStruct      Sampling structure
 nLoops            Number of iterations
 stepsPerPoint     Number of sampler steps per point saved
 maxtime           Amount of time to spend on calculation (in seconds)

OPTIONAL INPUTS
 proc              Number of processes if > 0.  Otherwise, the proces #.
 fdirectory        Do not use this parameter when calling function directly.  

OUTPUT
 sampleStruct      Sampling structure with sample points

 Jan Schellenberger 1/29/07
 (vaguely) based on code by:
 Markus Herrgard, Gregory Hannum, Ines Thiele, Nathan Price 4/14/06

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [sampleStruct] = ACHRSamplerParallelGeneral(sampleStruct,nLoops,stepsPerPoint, maxtime, proc, fdirectory)
0002 % ACHRSamplerParallelGeneral Artificial Centering Hit-and-Run sampler with in place (memory) point
0003 % managmenet
0004 %
0005 % sampleStruct = ACHRSamplerParallelGeneral(sampleStruct,nLoops,stepsPerPoint)
0006 %
0007 %INPUTS
0008 % sampleStruct      Sampling structure
0009 % nLoops            Number of iterations
0010 % stepsPerPoint     Number of sampler steps per point saved
0011 % maxtime           Amount of time to spend on calculation (in seconds)
0012 %
0013 %OPTIONAL INPUTS
0014 % proc              Number of processes if > 0.  Otherwise, the proces #.
0015 % fdirectory        Do not use this parameter when calling function directly.
0016 %
0017 %OUTPUT
0018 % sampleStruct      Sampling structure with sample points
0019 %
0020 % Jan Schellenberger 1/29/07
0021 % (vaguely) based on code by:
0022 % Markus Herrgard, Gregory Hannum, Ines Thiele, Nathan Price 4/14/06
0023 
0024 warning off MATLAB:divideByZero;
0025 
0026 
0027 %proc == 0 means master
0028 %proc greater than 1 means slave.
0029 if nargin < 5 % not parallel at all.
0030     parallel = 0;
0031     proc = 0;
0032 elseif nargin >= 5 && proc == 1 % not parallel (explicit)
0033     parallel = 0;
0034     proc = 0;
0035 else % parallel.
0036     parallel = 1;
0037     if proc < 0 % an indication that this is a slave process
0038         proc = -proc;
0039     else % indicator that you are a master process
0040         numproc = proc;
0041         proc = 0;
0042         % clear all files that may exist
0043         delete('xxMasterfil*.mat');
0044         delete('xxRound*.mat');
0045         delete('xxRoundDonePrint*.mat');
0046         delete('xxRoundAck*.mat');
0047         delete('xxDoneRound*.mat');
0048         delete('xxDoneP*.mat');
0049         delete('xxGlobalDone*.mat*');
0050     end
0051 end
0052 
0053 % Minimum allowed distance to the closest constraint
0054 maxMinTol = 1e-10;
0055 % Ignore directions where u is really small
0056 uTol = 1e-10; 
0057 safetycheck = false; % checks the direction of u for fixed directions.
0058 
0059 totalStepCount = 0;
0060 t0 = clock;
0061 
0062 
0063 if proc == 0 % if master thread
0064     if( ~ isfield(sampleStruct, 'points'))
0065         points = sampleStruct.warmupPts; % start with warmup points
0066     else
0067         points = sampleStruct.points; % continue with points
0068     end
0069     offset = sampleStruct.internal.offset;
0070 
0071     [dimX,nPoints] = size(points);
0072 
0073     points = points - offset*ones(1, nPoints);
0074 
0075     ub = sampleStruct.internal.ubnew;
0076     lb = sampleStruct.internal.lbnew;
0077     A = sampleStruct.internal.Anew;
0078     C = sampleStruct.internal.Cnew;
0079     D = sampleStruct.internal.Dnew;
0080     fixed = union(sampleStruct.internal.fixed,find(ub==lb));
0081     if (~isfield(sampleStruct.internal,'N'))
0082         if size(A,1)==0
0083             N=[];
0084             sampleStruct.internal.N=N;
0085         else
0086             if issparse(A)
0087                 N = null(full(A));
0088             else
0089                 N = null(A);
0090             end
0091             sampleStruct.internal.N = N;
0092         end
0093     else
0094         N = sampleStruct.internal.N;
0095     end
0096 
0097     movable = (1:dimX)';
0098     movable(fixed) = [];
0099     if safetycheck
0100         Nsmall = null(full(A(:, movable)));
0101     else
0102         Nsmall = [];
0103     end
0104 
0105     % Find the center of the space
0106     centerPoint = mean(points, 2);
0107 
0108     fidErr = fopen('ACHRParallelError.txt','w');
0109 
0110     pointRange = 1:nPoints;
0111     totalloops = nLoops;
0112 
0113     if parallel
0114         blah = 1;
0115         display('saving master file.');
0116         save('xxMasterfile', 'ub', 'lb', 'A', 'C', 'D', 'fixed', 'N', 'movable', 'Nsmall', 'numproc', 'nPoints' ); 
0117         display('finished saving master file.  spawning processes');
0118         for i = 1:(numproc - 1) % goes from 1 to 7 if proc == 8
0119             command = strcat('matlab -singleCompThread -automation -nojvm -r ACHRSamplerParallelGeneral([],',num2str(nLoops),',',num2str(stepsPerPoint),',0,', num2str(-i) ,',''' ,pwd, ''');exit; &' );
0120             display(command)
0121             system(command);
0122         end
0123         display('finished spawning processes');
0124     end
0125 end
0126 
0127 if proc > 0 %slave threads only
0128     cd (fdirectory);
0129     load('xxMasterfile', 'ub', 'lb', 'A', 'C', 'D', 'fixed', 'N', 'movable', 'Nsmall', 'numproc', 'nPoints');
0130     blah = 1;
0131 end
0132 
0133 for i = 1:nLoops
0134     if parallel % this whole block only gets executed in parallel mode.
0135         if proc == 0 %master thread does this.
0136             % save points
0137             display(strcat('distributing points round ', num2str(i)));
0138             save(strcat('xxRound', num2str(i)), 'points', 'centerPoint');
0139             save(strcat('xxRoundDonePrint', num2str(i)), 'blah');
0140             display(strcat('finished distributing points round ', num2str(i)));
0141         else  % if slave threads do this.
0142             display(strcat('reading in points round ', num2str(i)));
0143             while exist(strcat('xxRoundDonePrint', num2str(i), '.mat'), 'file') ~= 2; % wait for other thread to finish.
0144                 %display(strcat('waiting for round ', num2str(i)));
0145                 fprintf(1, '.');
0146                 if exist('xxGlobalDone.mat', 'file') == 2
0147                     exit;
0148                 end
0149                 pause(.25);
0150             end
0151             fprintf(1,'\nloading files four next round.\n');
0152             try load(strcat('xxRound', num2str(i)), 'points', 'centerPoint'); % load actual points
0153             catch pause(15) % for some reason at round 64 it needs extra time to load
0154                 load(strcat('xxRound', num2str(i)), 'points', 'centerPoint'); % load actual points
0155             end
0156             save(strcat('xxRoundAck', num2str(i),'x', num2str(proc) ), 'blah'); % save acknowledgement
0157             display(strcat('finished reading input and acknowledgment sent ', num2str(i)));
0158         end
0159         % divide up points.  master thread (proc = 0) gets first chunk.
0160         pointRange = subparts(nPoints, numproc, proc);
0161     end
0162     
0163     % actual sampling over pointRange
0164     for pointCount = pointRange
0165         % Create the random step size vector
0166         randVector = rand(stepsPerPoint,1);
0167         prevPoint = points(:,pointCount);
0168         curPoint = prevPoint;
0169         if mod(pointCount,200) == 0
0170             display(pointCount);
0171         end
0172         saveCoords = prevPoint(fixed);
0173         for stepCount = 1:stepsPerPoint
0174             % Pick a random warmup point
0175             randPointID = ceil(nPoints*rand);
0176             randPoint = points(:,randPointID);
0177             
0178             % Get a direction from the center point to the warmup point
0179             u = (randPoint-centerPoint);
0180             if ~isempty(fixed) % no need to reproject if there are no fixed reactions.
0181                 %ubefore = u;
0182                 if safetycheck
0183                     u(movable) = Nsmall * (Nsmall' * u(movable));
0184                 end
0185                 %uafter = u;
0186 
0187                 u(fixed) = 0; % takes care of biasing.
0188             end
0189             u = u/norm(u);
0190             
0191             % Figure out the distances to upper and lower bounds
0192             distUb = (ub - prevPoint);
0193             distLb = (prevPoint - lb);
0194             distD = (D-C*prevPoint);
0195             
0196             % Figure out positive and negative directions
0197             posDirn = (u > uTol);
0198             negDirn = (u < -uTol);
0199             move = C*u;
0200             posDirn2 = (move > uTol);
0201             negDirn2 = (move < -uTol);
0202             
0203             % Figure out all the possible maximum and minimum step sizes
0204             maxStepTemp = distUb./u;
0205             minStepTemp = -distLb./u;
0206             StepD = distD./move;
0207             maxStepVec = [maxStepTemp(posDirn);minStepTemp(negDirn);StepD(posDirn2 )];
0208             minStepVec = [minStepTemp(posDirn);maxStepTemp(negDirn);StepD(negDirn2 )];
0209         
0210             % Figure out the true max & min step sizes
0211             maxStep = min(maxStepVec);
0212             minStep = max(minStepVec);
0213             
0214             % Find new direction if we're getting too close to a constraint
0215             if (abs(minStep) < maxMinTol && abs(maxStep) < maxMinTol) || (minStep > maxStep)
0216                 fprintf('Warning small step: %f %f\n',minStep,maxStep);
0217                 continue;
0218             end
0219             
0220             % Pick a rand out of list_of_rands and use it to get a random
0221             % step distance
0222             stepDist = minStep + randVector(stepCount)*(maxStep-minStep);
0223             
0224             %fprintf('%d %d %d %f %f\n',i,pointCount,stepCount,minStep,maxStep);
0225             % Advance to the next point
0226             curPoint = prevPoint + stepDist*u;
0227             
0228             % Reproject the current point into the null space
0229             if mod (stepCount, 25) == 0
0230                 if ~isempty(N)
0231                     curPoint = N* (N' * curPoint);
0232                 end
0233                 curPoint(fixed) = saveCoords;
0234             end
0235 
0236             % Print out amount of constraint violation
0237             if (mod(totalStepCount,1000)==0) && proc == 0 % only do for master thread
0238               fprintf(fidErr,'%10.8f\t%10.8f\t',max(curPoint-ub),max(lb-curPoint));
0239             end
0240             
0241             % Move points inside the space if reprojection causes problems
0242             overInd = (curPoint > ub);
0243             underInd = (curPoint < lb);
0244             if (sum(overInd)>0) || (sum(underInd)>0)
0245               curPoint(overInd) = ub(overInd);
0246               curPoint(underInd) = lb(underInd);
0247             end
0248             
0249             % Print out amount of constraint violation
0250             if (mod(totalStepCount,1000) == 0) && proc == 0 % only do for master thread
0251               fprintf(fidErr,'%10.8f\n',full(max(max(abs(A*curPoint)))));
0252             end
0253             
0254             prevPoint = curPoint;
0255             
0256             % Count the total number of steps
0257             totalStepCount = totalStepCount + 1;
0258             
0259         end % Steps per point
0260         
0261         % Final reprojection
0262         if ~isempty(N)
0263             curPoint = N* (N' * curPoint);
0264         end
0265         curPoint(fixed) = saveCoords;
0266         centerPoint = centerPoint + (curPoint - points(:,pointCount))/nPoints; % only swapping one point... it's trivial.
0267         
0268         % Swap current point in set of points.
0269         points(:,pointCount) = curPoint;               
0270     end % Points per cycle
0271     
0272     if parallel % do this block if in parallel mode (regather points)
0273         if proc == 0 % if master
0274             % look for acknowledgements.
0275             display(strcat ('waiting for acknowledgement', num2str(i)));
0276             donewaiting = 0;
0277             while ~donewaiting
0278                 donewaiting = 1;
0279                 for k = 1:(numproc-1);
0280                     if exist(strcat('xxRoundAck',num2str(i), 'x', num2str(k),'.mat'), 'file') ~= 2
0281                         donewaiting = 0;
0282                     end
0283                 end
0284                 %display('waiting for acknowledgement');
0285                 fprintf(1, '.');
0286                 pause(.25)
0287             end
0288             % all other processes have received their information.  delete temporary files.
0289             fprintf(1, '\n');
0290             for k = 1:(numproc-1)
0291                 delete (strcat('xxRoundAck', num2str(i), 'x', num2str(k), '.mat'));
0292             end
0293             delete(strcat('xxRound', num2str(i), '.mat'));
0294             delete(strcat('xxRoundDonePrint', num2str(i), '.mat'));
0295             
0296             % look for return values.
0297             display(strcat ('waiting for return files ', num2str(i)));
0298             donewaiting = 0;
0299             while ~donewaiting
0300                 donewaiting = 1;
0301 
0302                 for k = 1:(numproc-1);
0303                     if exist(strcat('xxDoneP',num2str(i), 'x', num2str(k),'.mat'), 'file') ~= 2
0304                         donewaiting = 0;
0305                     end
0306                 end
0307                 pause(.25)
0308                 fprintf(1,'.');            
0309             end
0310             fprintf(1, '\nAll processes finished.  reading\n');
0311             for k = 1:(numproc-1)
0312                 load (strcat('xxDoneRound', num2str(i), 'x', num2str(k)), 'points2');
0313                 r2 = subparts(nPoints, numproc, k);
0314                 points(:,r2) = points2;
0315                 delete (strcat('xxDoneRound', num2str(i), 'x', num2str(k), '.mat') );
0316                 delete (strcat('xxDoneP',num2str(i), 'x', num2str(k), '.mat'));
0317             end
0318             centerPoint = mean(points, 2); % recalculate center point after gathering all data.
0319             display(strcat('done with round ', num2str(i)));
0320         else % if slave
0321             points2 = points(:, pointRange);
0322             save (strcat('xxDoneRound', num2str(i), 'x', num2str(proc)), 'points2');
0323             save (strcat('xxDoneP',num2str(i), 'x', num2str(proc)), 'blah');
0324         end
0325     end
0326     
0327     t1 = clock();
0328     fprintf('%10.0f s %d steps\n',etime(t1, t0),i*stepsPerPoint);
0329     if etime(t1, t0) > maxtime && proc == 0 % only master thread can terminate due to time limits.
0330         totalloops = i;
0331         break;
0332     end
0333 end
0334 
0335 if proc > 0 % slave threads terminate here.
0336     return;
0337 end
0338 points = points + offset*ones(1, nPoints);
0339 sampleStruct.points = points;
0340 if ~ isfield(sampleStruct, 'steps')
0341     sampleStruct.steps = 0;
0342 end
0343 
0344 sampleStruct.steps = sampleStruct.steps + stepsPerPoint*totalloops;
0345 % flag for all other handles to terminate.
0346 if parallel
0347     save('xxGlobalDone', 'blah');
0348     delete('xxMaster*.mat');
0349 end
0350 
0351 fclose(fidErr);
0352 
0353 function out = subparts(nPoints, n, k)
0354     out = (floor(nPoints*k/n)+1) : (floor(nPoints*(k+1)/n));
0355 return

Generated on Thu 21-Jun-2012 15:39:23 by m2html © 2003