Message Boards

GROUPS:

Staff Picks Computer Science Mathematica Add-Ons External Programs and Systems Graphics and Visualization Wolfram Language Packages Challenges Events & Media Quantum Computing

Nikolay Murzin

Solving Classiq Coding Competition with the Wolfram Quantum Framework

Nikolay Murzin

Posted 13 days ago

834 Views

4 Replies

17 Total Likes

Follow this post

enter image description here

Solving Classiq Coding Competition with the Wolfram Quantum Framework
by Nikolay Murzin

Recently tried myself at my first quantum competition and ended up winning one of the Innovative Solution prizes. Here is a slightly messy notebook in a form of a Computational Essay that formed my final submission:

Intro

The problem:

https://www.classiq.io/competition/kakuro

Kakuro is a logic puzzle, often referred to as a mathematical transliteration of the crossword. The puzzle is played on a grid of filled and empty cells. The filled cells contain the required sum for the matching empty cells.

The objective of the puzzle is to fill the empty spaces under two constraints:
- Two cells on the same row/column cannot have the same number.
- The sum of the cells on each row/column should equal the matching filled cell.

Using the QuantumFramework paclet, which is very much in early stages of development. Some features were written specifically to help solve this problem, like boolean oracles and improved Qiskit integration.

https://www.wolframcloud.com/obj/wolframquantumframework/DeployedResources/Paclet/Wolfram/QuantumFramework/

In[]:=

PacletInstall["https://wolfr.am/DevWQCF",ForceVersionInstall->True]<<Wolfram`QuantumFramework`

Out[]=

PacletObject

Name: Wolfram/QuantumFramework

Version: 1.0.16



In[]:=

mat=

x0	x1	0
x2	x3	x4
0	x5	x6

;

Manual solution for practice:

In[]:=

sol={x0->2,x1->1,x2->3,x3->2,x4->0,x5->0,x6->1}

Out[]=

{x02,x11,x23,x32,x40,x50,x61}

In[]:=

sol[[All,2]]

Out[]=

{2,1,3,2,0,0,1}

In[]:=

mat/.sol//MatrixForm

Out[]//MatrixForm=

2	1	0
3	2	0
0	0	1

Construct constraints:

In[]:=

constraints=Join[Catenate[Unequal@@@Subsets[DeleteCases[#,0],{2}]&/@mat],Catenate[Unequal@@@Subsets[DeleteCases[#,0],{2}]&/@(mat)],Thread[Plus@@@mat=={3,5,1}],Thread[Plus@@@(mat)=={5,3,1}]]

Out[]=

{x0≠x1,x2≠x3,x2≠x4,x3≠x4,x5≠x6,x0≠x2,x1≠x3,x1≠x5,x3≠x5,x4≠x6,x0+x13,x2+x3+x45,x5+x61,x0+x25,x1+x3+x53,x4+x61}

Allowed simplification

Constraints from the forum here and here:

Out[]=

(x0 != x1) and

(x0 != x2) and

(x1 != x3) and

(x1 != x5) and

(x2 + 2 != x3) and

(x2 + x4 + x3 == 3)

(x3 == 2) and

(x3 != x4) and

(x3 != x5) and

(x4 != x6) and

(x5 != x6) and

x3 = {x3[2], x3[1]}

In[]:=

boolSol=MapAt[ResourceFunction["InverseBoole"],{All,2}]@{x0->0,x1->1,x2->1,x3[1]->0,x3[2]->1,x4->0,x5->0,x6->1}

Out[]=

{x0False,x1True,x2True,x3[1]False,x3[2]True,x4False,x5False,x6True}

Semi-manual constraint construction:

Enumerate cases for x2 + x4 + x3 == 3

In[]:=

Cases[Tuples[{{0,1},{0,1},{0,1,2,3}}],tuple_/;Total[tuple]==3:>ResourceFunction["InverseBoole"]@MapAt[Splice@IntegerDigits[#,2,2]&,tuple,-1]]plusConstraint=Or@@(And@@Thread[Pick[{x2,x4,x3[2],x3[1]},#]]&/@%)

Out[]=

{{False,False,True,True},{False,True,True,False},{True,False,True,False},{True,True,False,True}}

Out[]=

(x3[2]&&x3[1])||(x4&&x3[2])||(x2&&x3[2])||(x2&&x4&&x3[1])

In[]:=

vars={x0,x1,x2,x3[1],x3[2],x4,x5,x6};boolConstraints={(x3[2]&&x3[1])||(x4&&x3[2])||(x2&&x3[2])||(x2&&x4&&x3[1])(*x2+x4+x3==3*),!x2&&(!x3[2]||x3[1])||x2&&(!x3[2]||!x3[1])(*x2+2!=x3*),!x4&&(x3[2]||x3[1])||x4&&(x3[2]||!x3[1])(*x3!=x4*),!x5&&(x3[2]||x3[1])||x5&&(x3[2]||!x3[1])(*x3!=x5*),!x3[1]&&x3[2](*x3==2*),x0!=x1,(*x0!=x1*)x0!=x2,(*x0!=x2*)x1!=x3[1],(*x1!=x3*)x1!=x5,(*x1!=x5*)x4!=x6(*x4!=x6*),x5!=x6(*x5!=x6*)}/.{Unequal[x_,y_]:>BooleanConvert[Not[x⧦y]],Equal[x_,y_]:>BooleanConvert[x⧦y]};

In[]:=

constraintLabels={"x2 + x4 + x3 == 3","x2 + 2 != x3","x3 != x4","x3 != x5","x3 == 2","x0 != x1","x0 != x2","x1 != x3","x1 != x5","x4 != x6","x5 != x6"};

In[]:=

And@@boolConstraints/.boolSol

Out[]=

True

Enumerate all, only single solution should satisfy all the constraints:

In[]:=

If[And@@boolConstraints/.Thread[vars->#],Thread[vars->Boole@#],Nothing]&/@Tuples[{False,True},Length[vars]]

Out[]=

Classical solver:

In[]:=

SatisfiabilityInstances[And@@boolConstraints,vars,All]

Out[]=

Circuit construction

Construct BooleanOracles for each individual constraint, putting each result on a separate ancilla:

In[]:=

oracle0=QuantumCircuitOperator[MapIndexed[QuantumCircuitOperator[{"BooleanOracle",#,vars,Length[vars]+First[#2]},constraintLabels[[First@#2]]]&,boolConstraints]];

In[]:=

oracle0["Flatten"]["Diagram"]

Out[]=

BooleanOracleR decomposes RZ(π) or RY(π) in place of an X with Gray Code Decomposition https://arxiv.org/pdf/quant-ph/0404089.pdf

In[]:=

oracleR=QuantumCircuitOperator[MapIndexed[QuantumCircuitOperator[{"BooleanOracleR",#,vars,Length[vars]+First[#2]},constraintLabels[[First@#2]]]&,boolConstraints]];

In[]:=

oracleR["Diagram",FontSize->8,"WireLabels"->None,ImageSize->Large,AspectRatio->1/2]

Out[]=

In[]:=

oracleR["Flatten"]["Diagram",FontSize->4,"WireLabels"->None,PlotRangePadding->0,ImageSize->Large,AspectRatio->1/16]

Out[]=

Counting gates by arity:

In[]:=

Counts[#["Arity"]&/@oracleR["Flatten"]["Operators"]]

Out[]=

151,284

84 CNOTs.

Define V-chain from Qiskit:

In[]:=

MCXCircuit // ClearAllMCXCircuit[order_List] := MCXCircuit[order] = Enclose @ Block[{controls = Most[order], target = Last[order], n = Length[order], max = Max[order], min = Min[order], circuit}, circuit = ExternalEvaluate[Wolfram`QuantumFramework`PackageScope`$PythonSession, "from qiskit import QuantumCircuit, transpilefrom qiskit.circuit.library import MCXVChain, MCXRecursivefrom wolframclient.language import wlimport picklenum_controls = <* n - 1 *>gate = MCXVChain(num_controls)num_ancillas = MCXVChain.get_num_ancilla_qubits(num_controls)qc = QuantumCircuit(<* max *> + num_ancillas)order = list(<* controls - 1 *>) + list(range(<* max *>, <* max *> + num_ancillas)) + [<* target - 1 *>]qc.append(gate, order, [])qc = transpile(qc, basis_gates=['cx', 'u1', 'u2', 'u3'], optimization_level=3)wl.Wolfram.QuantumFramework.QiskitCircuit(pickle.dumps(qc))"]; QuantumCircuitOperator[ConfirmMatch[circuit, _QiskitCircuit]["QuantumCircuit"], "V-chain"]]

In[]:=

DecomposeMCX // ClearAllDecomposeMCX[qo_QuantumOperator] /; MatchQ[qo["Label"], "Controlled"["NOT", _, _]] := Enclose @ Module[ controls = qo["ControlOrder"], target = qo["TargetOrder"], xs = Table[QuantumOperator["X", {i}], {i, qo["Label"][[3]]}], mcx, ConfirmAssert[Length[target] == 1]; mcx = MCXCircuit[Join[controls, target]]; If[ Length[xs] > 0, QuantumCircuitOperator[xs] @* mcx @* QuantumCircuitOperator[xs], mcx ]]DecomposeMCX[qo_QuantumOperator] := qoDecomposeMCX[qco_QuantumCircuitOperator] := QuantumCircuitOperator[DecomposeMCX /@ qco["Operators"]]

Try decomposing every MCX gate into a V-chain for comparison:

In[]:=

oracleVChain=DecomposeMCX@oracle0;

In[]:=

Counts[#["Arity"]&/@oracleVChain["Flatten"]["Operators"]]

Out[]=

1354,2225

The BooleanOracleR is much better... 84 is the final CNOT count without MCX

Going to use V-chain for the final MCX over 11 constraints with 60 CNOTs”

In[]:=

MCXCircuit[Range[m+1]]["Diagram",FontSize->4]

Out[]=

In[]:=

Count[MCXCircuit[Range[m+1]]["Operators"],op_/;MatchQ[op["Label"],"Controlled"[___]]]

Out[]=

Grover

Smaller test first:

In[]:=

SatisfiabilityInstances[(a||!b)&&(b||!c)&&(a&&!b||!d)&&(!b||d)&&(a&&!c),{a,b,c,d},All]

Out[]=

{{True,False,False,True},{True,False,False,False}}

In[]:=

n=4;m=5;

Test Grover circuit with oracle over 4 variables, 5 constraints and 2 solutions:

In[]:=

constraintsCircuit=QuantumCircuitOperator[{"BooleanOracleR",(a||b),{a,b,c,d},n+1,{"YRotation",Pi}}]/*QuantumCircuitOperator[{"BooleanOracleR",(b||!c),{a,b,c,d},n+2,{"YRotation",Pi}}]/*QuantumCircuitOperator[{"BooleanOracleR",(a&&!b||!d),{a,b,c,d},n+3,{"YRotation",Pi}}]/*QuantumCircuitOperator[{"BooleanOracleR",(!b||d),{a,b,c,d},n+4,{"YRotation",Pi}}]/*QuantumCircuitOperator[{"BooleanOracleR",(a&&!c),{a,b,c,d},n+5,{"YRotation",Pi}}];

In[]:=

groverTest=constraintsCircuit/*QuantumCircuitOperator[Table[QuantumOperator["H",{i}],{i,n+m}]]/*QuantumOperator[{"Controlled","NOT",Range[n+1,n+m]},{n+m+1}]/*QuantumCircuitOperator[Table[QuantumOperator["X",{i}],{i,n}]]/*QuantumOperator[{"Controlled","NOT",Range[n]},{n+m+1}]/*QuantumCircuitOperator[Table[QuantumOperator["X",{i}],{i,n}]]/*QuantumCircuitOperator[Table[QuantumOperator["H",{i}],{i,n}]];

In[]:=

groverTest["Diagram",FontSize->4]

Out[]=

V-chain version:

In[]:=

groverTestVChain=(*QuantumOperator["Z",{n+m+1}]/**)constraintsCircuit/*QuantumCircuitOperator[Table[QuantumOperator["H",{i}],{i,n+m}]]/*MCXCircuit[Range[n+1,n+m+1]]/*(*constraintsCircuit["Dagger"]/**)(*QuantumCircuitOperator[Table[QuantumOperator["H",{i}],{i,n+m}]]/**)(*QuantumOperator[{"Controlled","NOT",Range[n+1,n+m]},{n+m+1}]/**)QuantumCircuitOperator[Table[QuantumOperator["X",{i}],{i,n}]]/*QuantumOperator[{"Controlled","NOT",Range[n]},{n+m+1}]/*(*MCXCircuit[Range[n+1,n+m+1]]["Dagger"]/**)QuantumCircuitOperator[Table[QuantumOperator["X",{i}],{i,n}]]/*QuantumCircuitOperator[Table[QuantumOperator["H",{i}],{i,n}]];

In[]:=

groverTestVChain["Diagram",FontSize->4]

Out[]=

In[]:=

initStateTest=N@QuantumOperator["Z",{n+m+1}]@QuantumState[{"UniformSuperposition",n+m+1}];

Initial state with 3 ancillas, which are not in superposition:

In[]:=

initStateTestVChain0=QuantumState[{"Register",n+m+1+3,0}];

In[]:=

initStateTestVChain=N@QuantumOperator["Z",{n+m+1}]@QuantumOperator["H",Range[n+m+1,n+m+1+3]]@QuantumState[{"UniformSuperposition",n+m+1+3}];

Solutions are nicely recognisable after 2 iterations:

In[]:=

Nest[EchoFunction[QuantumMeasurementOperator[Range[n]]]@groverTest@#&,initStateTest,5]

QuantumMeasurement

Target: {1,2,3,4}

Measurement Outcomes: 16



QuantumMeasurement

Target: {1,2,3,4}

Measurement Outcomes: 16



QuantumMeasurement

Target: {1,2,3,4}

Measurement Outcomes: 16



QuantumMeasurement

Target: {1,2,3,4}

Measurement Outcomes: 16



QuantumMeasurement

Target: {1,2,3,4}

Measurement Outcomes: 16



Out[]=

QuantumState

Pure state	Qudits: 10
Type: Vector	Dimension: 1024
Picture: Schrödinger



Qiskit simulation for V-chain version:

In[]:=

groverTestQiskitN[k_]:=QuantumMeasurementOperator[Range[n]][RightComposition@@Table[groverTestVChain,k]]["Qiskit"]["Decompose",3]

In[]:=

groverTestQiskitWithInitN[k_]:=QuantumMeasurementOperator[Range[n]][QuantumCircuitOperator[{Splice@Table[QuantumOperator["H",{i}],{i,n+m}],QuantumOperator["Z",{n+m+1}]}]/*(RightComposition@@Table[groverTestVChain,k])]["Qiskit"]["Decompose",3]

In[]:=

Export[FileNameJoin[{NotebookDirectory[],"kakuro_test.qasm"}],groverTestQiskitWithInitN[3]["QASM"],"String"];

In[]:=

Map[BarChart[KeySort@groverTestQiskitN[#][initStateTestVChain],ChartLabelsPlaced[Automatic,Tooltip]]&,Range[6]]

Out[]=





Map[BarChart[KeySort@groverTestQiskitWithInitN[#][initStateTestVChain0],ChartLabelsPlaced[Automatic,Tooltip]]&,Range[6]]

Out[]=





Big circuit now (8 variables, 11 constraints, 1 solution):

In[]:=

n=Length[vars]m=Length[boolConstraints]

Out[]=

In[]:=

SatisfiabilityInstances[And@@boolConstraints,vars,All]

Out[]=

First without decomposition with just BooleanOracle.

Only gates corresponding to constraints (using RY(π) instead of a usual NOT):

In[]:=

constraintOracle=QuantumCircuitOperator[MapIndexed[QuantumCircuitOperator[{"BooleanOracle",#,vars,n+First[#2],{"YRotation",Pi}},constraintLabels[[First@#2]]]&,boolConstraints]];

In[]:=

constraintOracle["Operators"]

Out[]=

Adding Hadamards and MCX:

In[]:=

oracle0=constraintOracle/*QuantumCircuitOperator[Append[Table[QuantumOperator["H",{i}],{i,n+1,n+m}],QuantumOperator[{"Controlled","NOT",Range[n+1,n+m]},{n+m+1}]]]//N

Out[]=

In[]:=

initState=QuantumOperator["Z",{n+m+1}]@N@QuantumState[{"UniformSuperposition",n+m+1}];

In[]:=

grover=QuantumCircuitOperator[{"GroverDiffusion",n,QuantumOperator["NOT",{n+m+1}]},"Diffusion"]

Out[]=

The following runs slowly...

In[]:=

oracleState=oracle0@initState;//AbsoluteTiming

Out[]=

{33.4281,Null}

Test oracle:

In[]:=

test=QuantumTensorProduct[QuantumOperator[{"I",2,8}],QuantumOperator[QuantumState["111111111111"]["Dagger"],Range[n+1,n+m+1]]];

In[]:=

test[oracleState]["Formula"]

Out[]=

0.0441942|01101001〉

Diffuse two times:

In[]:=

finalState=grover@oracleState;//AbsoluteTiming

Out[]=

{33.3021,Null}

In[]:=

oracleState2=oracle0@finalState;//AbsoluteTiming

Out[]=

{38.2344,Null}

In[]:=

finalState2=grover@oracleState2;//AbsoluteTiming

Out[]=

{33.6321,Null}

Some optimised way to get just probabilities without QuantumMeasurement:

In[]:=

With[{proba=finalState2["Probabilities"]},KeyMap[Row]@TakeLargest[ResourceFunction["KeyGroupBy"][AssociationThread[IntegerDigits[Catenate@proba["ExplicitPositions"]-1,2,n+m+1],proba["ExplicitValues"]],Take[#,n]&,Total],3]]

Out[]=

011010010.0341975,001000100.00382488,000101000.00382488

The solution has the highest probability as expected.

Decomposed oracle with BooleanOracleR.

In[]:=

oracle=QuantumCircuitOperator[MapIndexed[QuantumCircuitOperator[{"BooleanOracleR",#,vars,n+First[#2],{"YRotation",Pi}},constraintLabels[[First@#2]]]&,boolConstraints]]/*QuantumCircuitOperator[{Splice@Table[QuantumOperator["H",{i}],{i,n+1,n+m}],QuantumOperator[{"Controlled","NOT",Range[n+1,n+m]},{n+m+1}]}]//N;

Not running it on any state, just confirming gate equality to the undecomposed version:

In[]:=

Thread[oracle0["Operators"]==oracle["Operators"]]

Out[]=

{True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True,True}

Replacing last MCXGate with V-chain to make the final oracle:

In[]:=

oracle["Operators"][[-1]]==QuantumOperator[{"Controlled","NOT",Range[n+1,n+m]},{n