DW2_code/000755 000765 000024 00000000000 12332735767 013415 5ustar00fspedalistaff000000 000000 DW2_code/._SEBREMforQUBO.m000644 000765 000024 00000000122 12276256154 016151 0ustar00fspedalistaff000000 000000 Mac OS X  2 RMATFMATLDW2_code/SEBREMforQUBO.m000644 000765 000024 00000054657 12276256154 015762 0ustar00fspedalistaff000000 000000 function [ RelativeEntropy,BestSolutions,BestObjective,BestSolutionFound,BestObjectiveFound ] = SEBREMforQUBO( Qfull,beta,Niterations,EmbeddingFlag,Step,Display) % SEBREMforQUBO Summary of this function goes here % Qfull is the QUBO matrix of the problem we want to embed % Parameters MAX_N_SOLUTIONS = 40000; Nqubits = 512; Nvariables = size(Qfull,1); params.num_reads = 10000; % Number of reads per run params.auto_scale = false; %params.annealing_time = 100; absJmax = 1; abshmax = 2; RelativeEntropy = zeros(1,Niterations); BestObjective = zeros(1,Niterations); %BestSolution = zeros(Nvariables,Niterations); %%%%%%% Start of code %%%%%%% % Generate and initial embedding for the QUBO [QubitMapping,couplers,VariableInteraction,h_chimera,J_chimera,hfull,Jfull] = InitialEmbedding(Qfull,EmbeddingFlag); % Define the function that computes the Ising energy on spin configurations IsingEnergy = @(S) (S'*Jfull*S + S'*hfull); % Start the iterative procedure for iter=1:Niterations % Solve the Ising model with quantum processor [ answer] = IsingConnectSolve(0.5*h_chimera,0.5*J_chimera,params); NumberOfSamples = size(answer.solutions,2); % Extract a distribution from the samples [SpinSolutions,ProbabilityOfSamples,EnergiesOfSamples{iter}] = ExtractDistribution(answer,QubitMapping); % Compute objective function on samples Gvec = GvecComputation(SpinSolutions,IsingEnergy); % Compute the values of the normalized original function on the samples G{iter} = Gvec; % Extract best solutions [BestObjective(iter),BestSolutions{iter}] = FindBestSolution(SpinSolutions,Gvec,Qfull); % Compute the relative entropy RelativeEntropy(iter) = -ProbabilityOfSamples*log(ProbabilityOfSamples') + beta*(Gvec*ProbabilityOfSamples'); % Compute the gradient of the relative entropy [GradRE_h,GradRE_J] = REGradient(SpinSolutions,ProbabilityOfSamples,Gvec,beta,VariableInteraction); Grad(iter,:) = [GradRE_h GradRE_J]; if iter > 1 GradInnerProduct = Grad(iter,:)*Grad(iter-1,:)'/(norm(Grad(iter-1,:))*norm(Grad(iter,:))); else GradInnerProduct = 0; end % Update Ising model if GradInnerProduct < -0 Step = Step/1.2; elseif GradInnerProduct > 0.5 Step = Step*1.3; end [h_chimera,J_chimera] = UpdateIsingModel( h_chimera,J_chimera,GradRE_h,GradRE_J,Step,QubitMapping,couplers,BestObjective,BestSolutions,iter); h{iter} = h_chimera; J{iter} = J_chimera; % Start printing information if Display == 1 DisplayInformation(RelativeEntropy,iter,BestObjective,BestSolutions,beta,Niterations,GradInnerProduct,NumberOfSamples); % Save data save SEBREMdata BestObjective BestSolutions G EnergiesOfSamples Grad RelativeEntropy h J ; end end % Extract best solution [BestObjectiveFound,bestindex] = min(BestObjective); BestSolutionFound = BestSolutions{bestindex}(:,1); end %% function [ QubitMapping,couplers,VariableInteraction,h_chimera,J_chimera,hfull,Jfull] = InitialEmbedding( Qfull,EmbeddingFlag ) %INITIALEMBEDDING Provides the starting embedding for SEBREM % It takes as input the full matrix of a QUBO, converts it to upper % diagonal (it assumes it is symmetric and discards the lower triangle), % converts it to an Ising model (h and J), and normalizes it. Then uses % some embedding heuristic controlled by EmbeddingFlag: % EmbeddingFlag = 1: Greedy embedding % EmbeddingFlag = 2: Stochastic greedy embedding % EmbeddingFlag = 3: Randomized direct embedding % % As output if provides the mapping between variables and qubits % (QubitMapping), a matrix with all the pairs of qubits assigned to % variables that share a coupler (couplers), and the same information but % written in terms of the variables those qubits are assigned. It also % provides the embedding as h_chimera and J_chimera, and the normalized % ising model as hfull and Jfull % Parameters Nqubits = 512; absJmax = 1; abshmax = 2; StartingQubit = 217; % First qubit of one of the central cells % Load adjacency matrix and qubits used load VesuviusAdjacencyMatrix Adjacency_matrix QubitsUsed; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Start of code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Convert Qfull to upper triangular Qfull = triu(Qfull) + triu(Qfull',1); % Convert to Ising model to allow proper normalization [hfull, Jfull,~] = quboToIsing(Qfull); % Now Jfull is also upper triangular % Normalize NormFactor = max([(1/abshmax)*max(abs(hfull)) (1/absJmax)*max(abs(Jfull))]); hfull = hfull/NormFactor; Jfull = Jfull/NormFactor; % Construct initial embedding if EmbeddingFlag == 1 % Apply the greedy embedding algorithm to find a starting embedding [ QubitMapping,AdjMat ] = GreedyEmbedding( Jfull,Adjacency_matrix,StartingQubit ); elseif EmbeddingFlag == 2 % Apply stochastic version of greedy embedding [ QubitMapping,AdjMat ] = StochasticGreedyEmbedding( Jfull,Adjacency_matrix,StartingQubit ); else [ QubitMapping,AdjMat ] = RandomizedDirectEmbedding( Jfull,Adjacency_matrix,QubitsUsed ); end % Embedded J: AdjMat is only an adjacency matrix % J_chimera should be upper triangular (because Jfull is) J_chimera = zeros(Nqubits); h_chimera = zeros(Nqubits,1); J_chimera(QubitMapping(:,2)',QubitMapping(:,2)') = Jfull.*AdjMat; J_chimera = triu(J_chimera); % Keep J_chimera upper triangular h_chimera(QubitMapping(:,2)') = hfull; % Assign local fields % Check if there is a zero row/column in J_chimera for which % h_chimera is also zero, and add a very small local field if that is the % case ZeroDetect = sum(abs(full([h_chimera(QubitMapping(:,2)') (sum(Jfull.*AdjMat + (Jfull.*AdjMat)'))'])),2); if ~isempty(find(ZeroDetect == 0,1)) h_chimera(QubitMapping(find(ZeroDetect == 0),2)) = 0.001; end % Inverse QubitMapping InverseQubitMapping = zeros(Nqubits,1); InverseQubitMapping(QubitMapping(:,2)) = QubitMapping(:,1); % Find all the available couplers between the qubits in Qchimera Reduced_Adjacency_matrix = zeros(Nqubits); Reduced_Adjacency_matrix(QubitMapping(:,2)',QubitMapping(:,2)') = Adjacency_matrix(QubitMapping(:,2)',QubitMapping(:,2)'); MaxCouplers = 0; for i=1:Nqubits nonzero_couplers = find(Reduced_Adjacency_matrix(i,1:Nqubits) == 1); for j=1:length(nonzero_couplers) if i < nonzero_couplers(j) MaxCouplers = MaxCouplers+1; couplers(MaxCouplers,1:2) = [i nonzero_couplers(j)]; % matrix with all the couplers between the qubits in QubitMaping VariableInteraction(MaxCouplers,1:2) = [InverseQubitMapping(i) ... InverseQubitMapping(nonzero_couplers(j))]; % Same as couplers but in terms of the original variables end % These couplers respect the upper triangular index ordering (couplers(ct,1) < couplers(ct,2)) end end couplers = couplers(1:MaxCouplers,:); % Keep only the assigned ones VariableInteraction = VariableInteraction(1:MaxCouplers,:); % Keep only the assigned ones end %% function [SpinSolutions,ProbabilityOfSamples,EnergiesOfSamples] = ExtractDistribution(answer,QubitMapping) %EXTRACTDISTRIBUTION Summary of this function goes here % Detailed explanation goes here SpinSolutions = answer.solutions(QubitMapping(:,2),:); % Keep only the values of the qubits associated with problem variables ProbabilityOfSamples = answer.num_occurrences/sum(answer.num_occurrences); % Estimate probability of each sample by its frequency EnergiesOfSamples = answer.energies; end %% function [ QubitMapping,AdjMat ] = GreedyEmbedding( J,AdjacencyMatrix,StartingQubit ) %GREEDYEMBEDDING Generates an embedding that prioritizes keeping the %strongest connections of J. % Inputs: % % J : Interaction matrix to be approximated % AdjacencyMatrix : Adjacency matrix of the chip % StartingQubit : Initial qubit to start the embedding (picking one with % high connectivity and somwhere around the center of the % chip recommended) % % % Outputs: % % QuibtMapping : N x 2 matrix, first column is a list of variables of % the quadratic function, second column has the qubits assigned to those % variables % AdjMat : reduced adjacency matrix representing the connections between % the qubits chosen from the embedding J = J - diag(diag(J)); Nvars = size(J,1); Nqubits = 512; % Total number of qubits in Vesuvius (used for indexing) Qubit = StartingQubit; MissingQubits = [75 123 141 361 367 382 392 396 424]; % Missing qubits in Vesuvius %Remove Missing qubits from AdjacencyMatrix (in case they're still there) for i=1:length(MissingQubits) AdjacencyMatrix(MissingQubits(i),:) = 0*AdjacencyMatrix(MissingQubits(i),:); AdjacencyMatrix(:,MissingQubits(i)) = 0*AdjacencyMatrix(:,MissingQubits(i)); end % MissingQubits = []; [~,StrongCouplingIndex] = max(sum(abs(J))); % Find the index with the strongest coupling VariablesLeft = 1:Nvars; VariablesAssigned = []; QubitsAssigned = []; QubitsLeft = 1:Nqubits; QubitsLeft(MissingQubits) = []; VariablesAssigned = StrongCouplingIndex; VariablesLeft(StrongCouplingIndex) = []; % Remove the assigned variable from list of variables J(:,StrongCouplingIndex) = zeros(Nvars,1); QubitsAssigned = [QubitsAssigned Qubit]; QubitNeighbors = find(AdjacencyMatrix(Qubit,:) ~= 0); % Qubits adjacent to first assigned qubit for n=2:Nvars for i=1:length(VariablesLeft) for j=1:length(QubitNeighbors) Coupling = 0; for k=1:length(VariablesAssigned) Coupling = Coupling + abs(J(VariablesAssigned(k),VariablesLeft(i)))*... AdjacencyMatrix(QubitNeighbors(j),QubitsAssigned(k)); end CouplingStrength(i,j) = Coupling; end end [MaxCouplingStrength,NewQubitIndex] = max(max(CouplingStrength)); %%% if size(CouplingStrength,1)> 1 PossibleNewQubitIndices = find(max(CouplingStrength) == MaxCouplingStrength); % Find all indices that give maximum coupling else PossibleNewQubitIndices = find(CouplingStrength == MaxCouplingStrength); % For the case CS is a row vector end for nqi=1:length(PossibleNewQubitIndices) DistanceToLastAssignedQubit(nqi) = abs(QubitsAssigned(end)-QubitNeighbors(PossibleNewQubitIndices(nqi))); end [~,MinIndex] = min(DistanceToLastAssignedQubit); % Find the qubit closest to last assigned qubit clear DistanceToLastAssignedQubit; NewQubit = QubitNeighbors(PossibleNewQubitIndices(MinIndex)); %%% %NewQubit = QubitNeighbors(NewQubitIndex); [~,NewVarIndex] = max(CouplingStrength(:,NewQubitIndex)); NewVar = VariablesLeft(NewVarIndex); clear CouplingStrength; VariablesAssigned = [VariablesAssigned NewVar]; VariablesLeft(find(VariablesLeft == NewVar)) = []; QubitsAssigned = [QubitsAssigned NewQubit]; QubitNeighbors = []; for k=1:length(QubitsAssigned) SingleQubitNeighbors = find(AdjacencyMatrix(QubitsAssigned(k),:) ~= 0); QubitNeighbors = unique([QubitNeighbors SingleQubitNeighbors]); end QubitNeighbors = setdiff(QubitNeighbors,QubitsAssigned); end QubitMapping = [VariablesAssigned' QubitsAssigned']; QubitMapping = sortrows(QubitMapping,1); AdjMat = AdjacencyMatrix(QubitMapping(:,2)',QubitMapping(:,2)'); end %% function [ Gvec] = GvecComputation( SpinSolutions, G ) %GVECCOMPUTATION Summary of this function goes here % Detailed explanation goes here CellSpinSolutions = mat2cell(SpinSolutions,size(SpinSolutions,1),ones(1,size(SpinSolutions,2))); G_CellSpinSolutions = cellfun(G,CellSpinSolutions,'UniformOutput',false); Gvec = cell2mat(G_CellSpinSolutions); end %% function [ distance ] = HammingDistance( String1,String2 ) %HAMMINGDISTANCE Summary of this function goes here % Detailed explanation goes here distance = sum(mod(String1+String2,2)); end %% function [GradRE_h,GradRE_J] = REGradient(SpinSolutions,ProbabilityOfSamples,Gvec,beta,VariableInteraction) %RECOMPUTATION COmputes the relative entropy and its gradient % It takes the structure answer as input and computes the relative % entropy between the input QUBO and the embedded one. Note that internally % the code works in the Ising framework, so solutions are strings of % {+1,-1}. Nvariables = size(SpinSolutions,1); % Fast gradient computation Corr1 = SpinSolutions*ProbabilityOfSamples'; % Vector of mean values of variables Corr2 = SpinSolutions*diag(ProbabilityOfSamples)*SpinSolutions'; VecCorr2 = vec(Corr2); % Vectorized two-variable correlation matrix Corr2coupler = VecCorr2(VariableInteraction(:,1) +(VariableInteraction(:,2) -1)*Nvariables); % Two variable correlations associated with just the couplers Hvec = ProbabilityOfSamples.*(log2(ProbabilityOfSamples) + beta*Gvec); ProdVec = SpinSolutions(VariableInteraction(:,1),:).*SpinSolutions(VariableInteraction(:,2),:); AverageLog = (log2(ProbabilityOfSamples) + beta*Gvec)*ProbabilityOfSamples'; GradRE_h = -beta*(SpinSolutions*Hvec' - Corr1*AverageLog)'; GradRE_J = -beta*(ProdVec*Hvec' - Corr2coupler*AverageLog)'; end %% function [h_chimera,J_chimera] = UpdateIsingModel( h_chimera,J_chimera,GradRE_h,GradRE_J,Step,QubitMapping,couplers,BestObjective,BestSolutions,iter) %UPDATEISINGMODEL Summary of this function goes here % Detailed explanation goes here absJmax = 1; abshmax = 2; Nvariables = size(QubitMapping,1); MaxCouplers = size(couplers,1); DELTA = Step/max(abs([GradRE_h GradRE_J])); for k=1:MaxCouplers %J_chimera(couplers(k,1),couplers(k,2)) = max(min(J_chimera(couplers(k,1),couplers(k,2)) - DELTA*GradRE_J(k) ,absJmax),-absJmax); % Add the random couplers: J_chimeraRandom stays upper triangular J_chimera(couplers(k,1),couplers(k,2)) = J_chimera(couplers(k,1),couplers(k,2)) - DELTA*GradRE_J(k) ; if k <= Nvariables %h_chimera(QubitMapping(k,2)) = max(min(h_chimera(QubitMapping(k,2)) - DELTA*GradRE_h(k),abshmax),-abshmax); % Add the random local fields h_chimera(QubitMapping(k,2)) = h_chimera(QubitMapping(k,2)) - DELTA*GradRE_h(k); end end Maxh = max(abs(h_chimera)); MaxJ = max(max(abs(J_chimera))); h_chimera = abshmax*h_chimera/Maxh; J_chimera = absJmax*J_chimera/MaxJ; MIXFACTOR = 0.0; MIXFACTOR = min(0.5,Step); if iter > 1 && BestObjective(iter) <= min(BestObjective(1:iter-1)) BinarySolutions = BestSolutions{iter}; [ h_mix,J_mix] = ChimeraParametersFromBinarySolution(BinarySolutions,QubitMapping,couplers); h_chimera = (1-MIXFACTOR)*h_chimera + max(abs(h_chimera))*MIXFACTOR*h_mix'; J_chimera = (1-MIXFACTOR)*J_chimera + max(max(abs(J_chimera)))*MIXFACTOR*J_mix; J_chimera = triu(J_chimera); % To keep J_chimera upper triangular end end %% function [ h,J] = ChimeraParametersFromBinarySolution(BinarySolutions,QubitMapping,couplers) %CHIMERAPARAMETERSFROMBINARYSOLUTION Summary of this function goes here % Detailed explanation goes here load VesuviusAdjacencyMatrix Adjacency_matrix QubitsUsed; Nqubits = size(QubitMapping,1); % Map BinarySolution to qubits and +-1 FullSolution = zeros(1,512); h = zeros(1,512); J = zeros(512); % Generate an ising model that combines all best solutions NumberOfBestSolutions = size(BinarySolutions,2); for soln=1:NumberOfBestSolutions FullSolution(QubitMapping(:,2)) = 2*BinarySolutions(:,soln)-1; h = h -(0.2/NumberOfBestSolutions)*FullSolution; % Check sign for i=1:size(couplers,1) J(couplers(i,1),couplers(i,2)) = J(couplers(i,1),couplers(i,2)) - (1/NumberOfBestSolutions)*FullSolution(couplers(i,1))*FullSolution(couplers(i,2)); end end end %% function [BestObjective,BestSolution] = FindBestSolution(SpinSolutions,Gvec,Qfull) % Detailed explanation goes here % Order Gvec [OrderedGvec,sortindices] = sortrows(Gvec'); % Order SpinSolutions OrderedSpinSolutions = SpinSolutions(:,sortindices); % Find how many solutions attain the lowest objective NumberOfGroundStates = length(find(OrderedGvec == OrderedGvec(1))); % Extract best solutions BestSolution = (OrderedSpinSolutions(:,1:NumberOfGroundStates) +1)/2; % Gives a binary string BestObjective = BestSolution(:,1)'*Qfull*BestSolution(:,1); end %% function DisplayInformation(RelativeEntropy,iter,BestObjective,BestSolutions,beta,Niterations,GradInnerProduct,NumberOfSamples) %DISPLAYINFORMATION Summary of this function goes here % Detailed explanation goes here if iter == 1 disp(' '); disp('SEQUENTIAL EMBEDDING BY RELATIVE ENTROPY MINIMIZATION '); disp(' '); fprintf('Parameters: beta = %f, \t Number of Iterations = %d \n', beta,Niterations); disp(' '); fprintf('Iter \t RE \t Best objective Hamming Distance \t Gradient angle \t Samples\n'); fprintf('%d \t %f \t %f \n',iter,RelativeEntropy(iter),BestObjective(iter)); else for i=1:size(BestSolutions{iter},2) for j=1:size(BestSolutions{iter-1},2) Hdistance(i,j) = HammingDistance(BestSolutions{iter-1}(:,j),BestSolutions{iter}(:,i)); end end hammingdistance = min(min(Hdistance)); fprintf('%d \t %f \t %f \t \t %d \t %f \t %d \n',iter,RelativeEntropy(iter),BestObjective(iter),hammingdistance,GradInnerProduct,NumberOfSamples); end end %% function [ QubitMapping,AdjMat ] = RandomizedDirectEmbedding( J,AdjacencyMatrix,QubitsUsed ) %RANDOMIZEDDIRECTEMBEDDING Summary of this function goes here % Detailed explanation goes here J = J - diag(diag(J)); Nvariables = size(J,1); Nqubits = 512; % Total number of qubits in Vesuvius (used for indexing) MissingQubits = [75 123 141 361 367 382 392 396 424]; % Missing qubits in Vesuvius %Remove Missing qubits from AdjacencyMatrix (in case they're still there) for i=1:length(MissingQubits) AdjacencyMatrix(MissingQubits(i),:) = 0*AdjacencyMatrix(MissingQubits(i),:); AdjacencyMatrix(:,MissingQubits(i)) = 0*AdjacencyMatrix(:,MissingQubits(i)); end % Create a random mapping form variables to qubits used QubitMapping(:,1) = 1:Nvariables; QubitMapping(:,2) = QubitsUsed(randperm(Nvariables))'; % Compute adjacency matrix of mapped variables AdjMat = AdjacencyMatrix(QubitMapping(:,2),QubitMapping(:,2)); end %% function [ QubitMapping,AdjMat ] = StochasticGreedyEmbedding( J,AdjacencyMatrix,StartingQubit ) %GREEDYEMBEDDING Generates an embedding that prioritizes keeping the %strongest connections of J. % Inputs: % % J : Interaction matrix to be approximated % AdjacencyMatrix : Adjacency matrix of the chip % StartingQubit : Initial qubit to start the embedding (picking one with % high connectivity and somwhere around the center of the % chip recommended) % % % Outputs: % % QuibtMapping : N x 2 matrix, first column is a list of variables of % the quadratic function, second column has the qubits assigned to those % variables % AdjMat : reduced adjacency matrix representing the connections between % the qubits chosen from the embedding J = J - diag(diag(J)); Nvars = size(J,1); Nqubits = 512; % Total number of qubits in Vesuvius (used for indexing) Qubit = StartingQubit; MissingQubits = [75 123 141 361 367 382 392 396 424]; % Missing qubits in Vesuvius %Remove Missing qubits from AdjacencyMatrix (in case they're still there) for i=1:length(MissingQubits) AdjacencyMatrix(MissingQubits(i),:) = 0*AdjacencyMatrix(MissingQubits(i),:); AdjacencyMatrix(:,MissingQubits(i)) = 0*AdjacencyMatrix(:,MissingQubits(i)); end % MissingQubits = []; [~,StrongCouplingIndex] = max(sum(abs(J))); % Find the index with the strongest coupling VariablesLeft = 1:Nvars; VariablesAssigned = []; QubitsAssigned = []; QubitsLeft = 1:Nqubits; QubitsLeft(MissingQubits) = []; VariablesAssigned = StrongCouplingIndex; VariablesLeft(StrongCouplingIndex) = []; % Remove the assigned variable from list of variables J(:,StrongCouplingIndex) = zeros(Nvars,1); QubitsAssigned = [QubitsAssigned Qubit]; QubitNeighbors = find(AdjacencyMatrix(Qubit,:) ~= 0); % Qubits adjacent to first assigned qubit for n=2:Nvars for i=1:length(VariablesLeft) for j=1:length(QubitNeighbors) Coupling = 0; for k=1:length(VariablesAssigned) Coupling = Coupling + abs(J(VariablesAssigned(k),VariablesLeft(i)))*... AdjacencyMatrix(QubitNeighbors(j),QubitsAssigned(k)); end CouplingStrength(i,j) = Coupling; end end [MaxCouplingStrength,NewQubitIndex] = max(max(CouplingStrength)); %%% if size(CouplingStrength,1)> 1 PossibleNewQubitIndices = find(max(CouplingStrength) == MaxCouplingStrength); % Find all indices that give maximum coupling else PossibleNewQubitIndices = find(CouplingStrength == MaxCouplingStrength); % For the case CS is a row vector end for nqi=1:length(PossibleNewQubitIndices) DistanceToLastAssignedQubit(nqi) = abs(QubitsAssigned(end)-QubitNeighbors(PossibleNewQubitIndices(nqi))); end [~,MinIndex] = min(DistanceToLastAssignedQubit); % Find the qubit closest to last assigned qubit clear DistanceToLastAssignedQubit; NewQubit = QubitNeighbors(PossibleNewQubitIndices(randi(length(PossibleNewQubitIndices)))); % Randomly choose next qubit %%% %NewQubit = QubitNeighbors(NewQubitIndex); [~,NewVarIndex] = max(CouplingStrength(:,NewQubitIndex)); NewVar = VariablesLeft(NewVarIndex); clear CouplingStrength; VariablesAssigned = [VariablesAssigned NewVar]; VariablesLeft(find(VariablesLeft == NewVar)) = []; QubitsAssigned = [QubitsAssigned NewQubit]; QubitNeighbors = []; for k=1:length(QubitsAssigned) SingleQubitNeighbors = find(AdjacencyMatrix(QubitsAssigned(k),:) ~= 0); QubitNeighbors = unique([QubitNeighbors SingleQubitNeighbors]); end QubitNeighbors = setdiff(QubitNeighbors,QubitsAssigned); end QubitMapping = [VariablesAssigned' QubitsAssigned']; QubitMapping = sortrows(QubitMapping,1); AdjMat = AdjacencyMatrix(QubitMapping(:,2)',QubitMapping(:,2)'); end %% function x = vec(X) [m n] = size(X); x = reshape(X,m*n,1); endDW2_code/._Solve0_1ILP.m000644 000765 000024 00000000122 12332464243 015724 0ustar00fspedalistaff000000 000000 Mac OS X  2 RMATFMATLDW2_code/Solve0_1ILP.m000644 000765 000024 00000002261 12332464243 015515 0ustar00fspedalistaff000000 000000 function [ ILPmatrix,VariableMap,ILPsolution,ILPobjective,SplitVariables,OutputString ] = Solve0_1ILP( filename,Penalty,beta,Niterations,EmbeddingFlag,Step,Display ) %SOLVE0_1ILP Solves 0-1 ILP for Minimal Separating Set using SEBREM % Filename contains the ILP description. It should be a .txt file, each % line a sequence of integers representing the variable indices present % in the equation described by that line % Read the filenmae and construct the ILP matrix [ ILPmatrix,VariableMap ] = ReadILP( filename ); % Convert the ILP matrix into a QUBO matrix adding slack variables A = ILPmatrix; f = ones(1,size(A,2)); b = ones(size(A,1),1); [ QUBO_matrix] = ILPtoQUBO( f,A,b,Penalty ); % Solve QUBO problem using SEBREM [ RelativeEntropy,BestSolutions,BestObjective,BestSolutionFound,BestObjectiveFound ] = SEBREMforQUBO( QUBO_matrix,beta,Niterations,EmbeddingFlag,Step,Display); % Extract solution for non-slack variables ILPsolution = BestSolutionFound(1:length(f)); ILPobjective = f*ILPsolution; VariableIndices = VariableMap(:,2)'; SplitVariables = VariableIndices(find(ILPsolution == 1)); % Format output OutputString = int2str(SplitVariables); end DW2_code/._SolverILP.m000644 000765 000024 00000000122 12332464201 015600 0ustar00fspedalistaff000000 000000 Mac OS X  2 RMATFMATLDW2_code/SolverILP.m000644 000765 000024 00000000675 12332464201 015400 0ustar00fspedalistaff000000 000000 function SolverILP( filename ) %SOLVERILP Takes an ILP in 'filename', solves it and prints split variables %to stdout Penalty = 10; beta=1; Niterations = 2; EmbeddingFlag = 3; Step = 1; Display = 0; [ ILPmatrix,VariableMap,ILPsolution,ILPobjective,SplitVariables,OutputString ] = Solve0_1ILP( filename,Penalty,beta,Niterations,EmbeddingFlag,Step,Display ); % Format output and print to stdout fprintf(1,OutputString); exit; end