Grape optimization

.gitignore

@@ -0,0 +1,4 @@
+# Eclipse project file
+
+.project
+

src/GRAPE.m

@@ -545,28 +545,15 @@
 (*Takes a pulse and finds the overall unitary.*)
 
 
-PropagatorFromPulse[pulse_,Hint_,Hcontrol_]:=
-	Module[{step=1,dt,dts,amps},
-		{dts,amps} = SplitPulse[pulse];
-		(* Notice the delay equal on dt! Not a bug. *)
-		dt := dts[[step++]];
-		Fold[
-			MatrixExp[-I(Hint+#2.Hcontrol) dt].#1&,
-			IdentityMatrix[Length[Hint]],
-			amps
-		]
-	]
-
-
 PropagatorListFromPulse[pulse_,Hint_,Hcontrol_]:=
-	Module[{dts,dt,step=1,amps},
-		{dts,amps} = SplitPulse[pulse];
-		dt := dts[[step++]];
-		Map[
-			MatrixExp[-I(Hint+#.Hcontrol)dt]&,
-			amps
-		]
-	]
+	Map[
+		MatrixExp[-I* (Hint+Rest[#].Hcontrol)* First[#]]&,
+		pulse
+	];
+
+PropagatorFromPulse[pulse_,Hint_,Hcontrol_]:= Fold[
+	Dot, PropagatorListFromPulse[pulse, Hint, Hcontrol]
+];
 
 
 (* ::Subsubsection::Closed:: *)
@@ -2517,23 +2504,22 @@
 
 				(********************** CHOOSE A DIRECTION ***********************)
 
-				Module[{distortedPulse, distortionJacobian},
+				Module[{distortedPulse, distortionJacobian, calculateUtilities},
 
 					(* Distort the current pulse and calculate the Jacobian with the control knobs. *)
 					(* The output may have distribution symbols included *)
+					Print["Calculating Initial Distortion"];
 					If[\[Not]distortionDependsOnDist,
 						{distortedPulse, distortionJacobian} = DistortionFn[pulse, True]
 					];
 
-					(* Now we do a big weighted sum over all elements of the distribution *)
-					{rawUtility, utility, gradient} = Sum[
-						Module[
+					(* Calculates the gradient for a given rule *)
+					calculateUtilities[parameterProbabilityAndRule_] := Module[
 							{reps, prob, objFunVal, objFunGrad, penaltyGrad, totalGrad},
 
 							(* The replacements and probability of the current distribution element *)
-							reps=distReps[[d]];
-							prob=distPs[[d]];
-
+							reps=parameterProbabilityAndRule[[2]];
+							prob=parameterProbabilityAndRule[[1]];
 							(* If the distortion depends on the distribution, we evaluate it here, replacing what we can.
 							Anything else will get replaced in the call to UtilityGradient below. *)
 							If[distortionDependsOnDist,
@@ -2553,10 +2539,15 @@
 
 							(* Update the utility with the penalty *)
 							prob*{objFunVal, objFunVal - penalty, totalGrad}
-						],
-						{d,distNum}
-					];
-
+						];
+					Print["Calculating Utilities to Determine Jacobian"];
+					{rawUtility, utility, gradient} = Total[
+						ParallelMap[
+							calculateUtilities, 
+							Transpose[{distPs, distReps}], 
+							Method -> "CoarsestGrained"]
+						];
+
 					(* Apply the mask to the gradient, so we can avoid bad controls. *)
 					gradient = derivMask * badControlLimitGradientMask * gradient;
 				];
@@ -2639,16 +2630,25 @@
 
 				(************** CALCULATE STEP SIZE MULTIPLIER **************)
 				(* Evaluate the objective function at two points along the gradient direction *)
-				Module[{testFn},
-					testFn = With[{distortedPulse=DistortionFn[pulse+AddTimeSteps[0,#*stepSize*(\[Epsilon]max^2*#&/@goodDirec)], False]},
-							Sum[
-								With[{reps=distReps[[d]], prob=distPs[[d]]},
-									prob*(Utility[PropagatorFromPulse[distortedPulse/.reps,Hint/.reps,Hcontrol/.reps], target/.reps] - 
-										PulsePenaltyFn[distortedPulse/.reps, False])
-								],
-								{d, distNum}
-							]
-						]&;
+				Module[{testFn, calculateUtilityForReplacement},
+
+					calculateUtilityForReplacement[probabilityAndRule_] := Function[{distortedPulse}, 
+						With[{reps=probabilityAndRule[[2]], prob=probabilityAndRule[[1]]},
+							prob*(Utility[PropagatorFromPulse[distortedPulse/.reps,Hint/.reps,Hcontrol/.reps], target/.reps] - 
+								PulsePenaltyFn[distortedPulse/.reps, False])
+						]];
+
+					Print["Calculating Distortion for Step Size"];
+
+					testFn = With[
+						{
+							distortedPulse=DistortionFn[pulse+AddTimeSteps[0,#*stepSize*((\[Epsilon]max^2*#)&/@goodDirec)], False]
+						},
+						Total[
+							ParallelMap[(calculateUtilityForReplacement[#][distortedPulse])&,
+							Transpose[{distPs, distReps}], 
+							Method -> "CoarsestGrained"]
+						]]&;
 					bestMult = lineSearchMeth[utility, testFn];
 				];
 
@@ -2660,6 +2660,8 @@
 
 				(* Make sure that the pulse does not exceed the limits *)
 				pulse=OptionValue[PulseLegalizer][pulse,controlRange];
+
+				Print["Updating Pulse"];
 
 				(* Pulse has changed; update best pulse for monitor *)
 				UpdateBestPulse;

test/GRAPETests.m

@@ -22,7 +22,7 @@
 (*THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THEIMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE AREDISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLEFOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIALDAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS ORSERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVERCAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USEOF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.*)
 
 
-(* ::Subsection::Closed:: *)
+(* ::Subsection:: *)
 (*Preamble*)
 
 
@@ -34,6 +34,7 @@
 
 
 Needs["QUDevTools`"];
+Needs["QuantumSystems`"];
 Needs["GRAPE`"];
 
 
@@ -51,15 +52,99 @@
 End[];
 
 
-(* ::Section::Closed:: *)
+(* ::Section:: *)
 (*Unit Tests*)
 
 
 Begin["`UnitTests`"];
 
-
-(* ::Subsection::Closed:: *)
-(*End*)
+(* ::Subsection:: *)
+(*PropagatorListFromPulse*)
+TestCase[$RegisteredTests, "PropagatorListFromPulse:AssertListCorrect",
+	Module[{randomPulse},
+		randomPulse = RandomPulse[1, 10, {{-1, 1}}];
+		Length[FromPulse[randomPulse]] == Length[PropagatorListFromPulse[randomPulse]]
+	]
+]
+
+(* ::Subsection:: *)
+(*Pulse Object*)
+
+
+(* ::Subsubsection:: *)
+(*Initialization*)
+
+
+TestCase[$RegisteredTests, "PulseDataStructure:InitializeHControl", 
+	Module[{p, HControl},
+		HControl = {RandomHermitian[8], RandomHermitian[8]};		
+		p = Pulse[
+			ControlHamiltonians -> HControl
+		];
+		ControlHamiltonians[p] == Hinternal;
+	]
+]
+TestCase[$RegisteredTests, "PulseDataStrucutre:InitializeHInternal",
+	Module[{p, Hinternal},
+	Hinternal = RandomHermitian[8];
+	p = Pulse[InternalHamiltonian -> Hinternal];
+	InternalHamiltonian[p] == Hinternal;
+]]
+
+TestCase[$RegisteredTests, "PulseDataStructure:InitializeUtarget",
+	Module[{p, target},
+		target = RandomUnitary[8];
+		p = Pulse[
+			Target -> target
+		];
+		Target[p] == target
+	]
+]
+
+
+(* ::Subsection:: *)
+(*FindPulse*)
+
+
+(* ::Subsubsection:: *)
+(*Simple Pulse*)
+
+
+TestCase[$RegisteredTests, "FindPulse:SimplePulse",
+	Module[{
+		pulse, Hint, HControl, controlRange, initialGuess, target, \[Phi]target
+	},
+		Hint = TP[Z];
+		HControl = 2\[Pi] * {TP[X], TP[Y]};
+		controlRange = {{-1, 1}, {-1, 1}};
+		target = RandomUnitary[2];
+		\[Phi]target = 0.99;
+		initialGuess = RandomSmoothPulse[1, 10, controlRange];
+		pulse = FindPulse[initialGuess, target, \[Phi]target, controlRange, HControl, Hint];
+		Re[Tr[pulse[Target] . target / 2]] >= \[Phi]target;
+]
+]
+
+
+(* ::Subsubsection:: *)
+(*Repetitions*)
+
+
+TestCase[$RegisteredTests, "FindPulse:Repetitions",
+	Module[{
+		pulse, repetitions, Hint, HControl, controlRange, initialGuess, target, \[Phi]target
+	},
+		Hint = TP[Z];
+		repetitions = 10;
+		HControl = 2\[Pi] * {TP[X], TP[Y]};
+		controlRange = {{-1, 1}, {-1, 1}};
+		target = RandomUnitary[2];
+		\[Phi]target = 0.99;
+		initialGuess = RandomSmoothPulse[1, 10, controlRange];
+		pulse = FindPulse[initialGuess, target, \[Phi]target, controlRange, HControl, Hint, Repetitions -> repetitions];
+		Re[Tr[pulse[Target] . target / 2]] >= \[Phi]target;
+]
+]
 
 
 End[];