{-# LANGUAGE NoMonomorphismRestriction #-} -- | This module contains an implementation of the oracle and its -- automatic lifting to quantum circuits using Template Haskell. module Quipper.Algorithms.QLS.TemplateOracle where import Data.Complex import Quipper.Utils.Auxiliary import Control.Monad import Quipper.Libraries.Arith hiding (template_symb_plus_) import Quipper.Algorithms.QLS.Utils import Quipper.Algorithms.QLS.QDouble import Quipper.Algorithms.QLS.RealFunc import Quipper.Algorithms.QLS.QSignedInt import Quipper.Algorithms.QLS.CircLiftingImport import Quipper import Quipper.Internal.CircLifting -- | Lifted version of @'any'@ (local version). build_circuit local_any :: (a -> Bool) -> [a] -> Bool local_any f l = case l of [] -> False (h:t) -> if (f h) then True else local_any f t -- | Auxiliary function. build_circuit itoxy :: Int -> Int -> Int -> (Int, Int) itoxy i nx ny = if (i <= 0) then (0,0) else if (i > (nx-1)*ny + nx*(ny-1)) then (0,0) else ((mod (i - 1) (2*nx - 1)) + 1, ceiling ((fromIntegral i)/(2.0*(fromIntegral nx)-1.0))) -- | The function sin /x/ \/ /x/. build_circuit sinc :: Double -> Double sinc x = if (x /= 0.0) then (sin x) / x else 1.0 -- | Auxiliary function. build_circuit edgetoxy :: Int -> Int -> Int -> (Double, Double) edgetoxy e nx ny = let (ex,ey) = itoxy e nx ny in if (ex < nx) then ((fromIntegral ex) + 0.5, fromIntegral ey) else (fromIntegral (ex - nx + 1), (fromIntegral ey) + 0.5) -- | Auxiliary function. The inputs are: -- -- * /y1/ :: 'Int' - Global edge index of row index of desired matrix element; -- -- * /y2/ :: 'Int' - Global edge index of column index of desired matrix element; -- -- * /nx/ :: 'Int' - Number of vertices left to right; -- -- * /ny/ :: 'Int' - Number of vertices top to bottom; -- -- * /lx/ :: 'Double' - Length of horizontal edges (distance between vertices in /x/ direction); -- -- * /ly/ :: 'Double' - Length of vertical edges (distance between vertices in /y/ direction); -- -- * /k/ :: 'Double' - Plane wave wavenumber. -- -- The output is the matrix element /A/(/y1/, /y2/). build_circuit calcmatrixelement :: -- -- Inputs Int -- int y1 - Global edge index of row index of desired matrix element -> Int -- int y2 - Global edge index of column index of desired matrix element -> Int -- int nx - Number of vertices left to right -> Int -- int ny - Number of vertices top to bottom -> Double -- float lx - Length of horizontal edges (distance between vertices in x direction) -> Double -- float ly - Length of vertical edges (distance between vertices in y direction) -> Double -- float k - Plane wave wavenumber -- Outputs -> Complex Double -- float A - Matrix element A(y1,y2) calcmatrixelement y1 y2 nx ny lx ly k = let (xg1, yg1) = itoxy y1 nx ny in let (xg2, yg2) = itoxy y2 nx ny in let b = if ( (y1==y2) && (xg1 >= nx) ) then ly/lx - k*k*lx*ly/3.0 -- B11 and B22 else if ( (y1==y2) && (xg1=nx) ) then -ly/lx + k*k*lx*ly/12.0 -- B13 else if ( (yg1==(yg2+1)) && (xg1==(xg2-nx+1)) && (xg2>=nx) ) then -1.0 -- B31 else if ( (yg2==(yg1+1)) && (xg2==(xg1-nx+1)) && (xg1>=nx) ) then -1.0 -- B14 else if ( (yg1==(yg2+1)) && (xg1==(xg2-nx)) && (xg1=nx) ) then -1.0 -- B24 else if ( (yg2==yg1) && (xg2==(xg1+nx)) && (xg2>=nx) ) then -1.0 -- B32 else if ( (yg1==yg2) && (xg1==(xg2+nx-1)) && (xg2 [] (x:tt) -> case tt of [] -> [] (y:t) -> (x,y):(get_edges (y:t)) -- | Auxiliary function. build_circuit checkedge :: Int -> [(Double,Double)] -> Int -> Int -> Bool checkedge e scatteringnodes nx ny = let (xi,yi) = edgetoxy e nx ny in let test_elt ((x1,y1),(x2, y2)) = xi >= x1 && yi >= y1 && xi <= x2 && yi <= y2 in let half = take_half scatteringnodes in let test_list = local_any test_elt (get_edges half) in (not test_list) -- | Oracle /r/. build_circuit calcRweights :: Int -> Int -> Int -> Double -> Double -> Double -> Double -> Double -> Complex Double calcRweights y nx ny lx ly k theta phi = let (xc',yc') = edgetoxy y nx ny in let xc = (xc'-1.0)*lx - ((fromIntegral nx)-1.0)*lx/2.0 in let yc = (yc'-1.0)*ly - ((fromIntegral ny)-1.0)*ly/2.0 in let (xg,yg) = itoxy y nx ny in if (xg == nx) then let i = (mkPolar ly (k*xc*(cos phi)))* (mkPolar 1.0 (k*yc*(sin phi)))* ((sinc (k*ly*(sin phi)/2.0)) :+ 0.0) in let r = ( cos(phi) :+ k*lx )*((cos (theta - phi))/lx :+ 0.0) in i * r else if (xg==2*nx-1) then let i = (mkPolar ly (k*xc*cos(phi)))* (mkPolar 1.0 (k*yc*sin(phi)))* ((sinc (k*ly*sin(phi)/2.0)) :+ 0.0) in let r = ( cos(phi) :+ (- k*lx))*((cos (theta - phi))/lx :+ 0.0) in i * r else if ( (yg==1) && (xg Int -> Int -> Int -> Int convertband y b nx ny = let nedges = (nx - 1)*ny + nx*(ny - 1) in let (ex,ey) = itoxy y nx ny in let x = if ( (ex < nx) && (ey /= 1) ) then case b of 1 -> y-2*nx+1 2 -> y-nx 3 -> y-nx+1 5 -> y 7 -> y+nx-1 8 -> y+nx 9 -> y+2*nx-1 _ -> -1 else if ( (ex < nx) && (ey == 1) ) then case b of 5 -> y 7 -> y+nx-1 8 -> y+nx 9 -> y+2*nx-1 _ -> -1 else if ( (ex >= nx) && (ex /= nx) && (ex /= 2*nx-1) ) then case b of 2 -> y-nx 3 -> y-nx+1 4 -> y-1 5 -> y 6 -> y+1 7 -> y+nx-1 8 -> y+nx _ -> -1 else if ( (ex >= nx) && (ex == nx) ) then case b of 3 -> y-nx+1 5 -> y 6 -> y+1 8 -> y+nx _ -> -1 else if ( (ex >= nx) && (ex == 2*nx-1) ) then case b of 2 -> y-nx 4 -> y-1 5 -> y 7 -> y+nx-1 _ -> -1 else -1 in if ( (x < 1) || (x > nedges) ) then -1 else x -- | Oracle /A/. It is equivalent to the Matlab function -- /getBandNodeValues/. -- -- 'getNodeValuesMoreOutputs' /v/ /b/ ... outputs the node of the -- edge connected to vertex /v/ in band /b/, and a real number -- parameterized by the 'BoolParam' parameter: the magnitude -- ('PFalse') or the phase ('PTrue') of the complex value at the -- corresponding place in the matrix /A/. build_circuit getNodeValuesMoreOutputs :: Int -> Int -> Int -> Int -> [(Double,Double)] -> Double -> Double -> Double -> BoolParam -> Int -> (Int, Double) getNodeValuesMoreOutputs v' b' nx ny scatteringnodes lx ly k argflag maxConnectivity = let maxC = getIntFromParam maxConnectivity in let b = getIntFromParam b' in let nedges = (nx - 1)*ny + nx*(ny - 1) in let flag = v' <= nedges in let v = (if flag then v' else (v' - nedges)) in let nodeDefault = (if flag then v + nedges else v) in let valueDefault = 0.0 in let indicesDefault = (-1, -1) in let isvalid = checkedge v scatteringnodes nx ny in if ( (not isvalid) && b == 5 ) then let indices = (if flag then (v, v + nedges) else (v + nedges, v)) in case argflag of PTrue -> (nodeDefault, valueDefault) PFalse -> (nodeDefault, 1.0) else if ( (not isvalid) && b /= 5) then (nodeDefault, valueDefault) else if ((b > maxC + 2) || (b <= 0)) then (nodeDefault, valueDefault) else let x = convertband v b' nx ny in if (x == -1) then (nodeDefault, valueDefault) else let isvalid = checkedge x scatteringnodes nx ny in if isvalid then let ax = calcmatrixelement v x nx ny lx ly k in let (node, indices) = if flag then (x + nedges, (v, x + nedges)) else (x, (v+nedges,x)) in let value = case argflag of PTrue -> atan2 (imagPart ax) (realPart ax) PFalse -> magnitude ax in (node, value) else (nodeDefault, valueDefault) -- | Auxiliary function for oracle /b/. The inputs are: -- -- * /y/ :: 'Int' - Global edge index. Note this is the unmarked /y/ -- coordinate, i.e. the coordinate without scattering regions removed; -- -- * /nx/ :: 'Int' - Number of vertices left to right; -- -- * /ny/ :: 'Int' - Number of vertices top to bottom; -- -- * /lx/ :: 'Double' - Length of horizontal edges (distance between vertices in /x/ direction); -- -- * /ly/ :: 'Double' - Length of vertical edges (distance between vertices in /y/ direction); -- -- * /k/ :: 'Double' - Plane wave wavenumber; -- -- * θ :: 'Double' - Direction of wave propagation; -- -- * /E0/ :: 'Double' - Magnitude of incident plane wave. -- -- The output is the magnitude of the electric field on edge /y/. build_circuit calcincidentfield :: -- -- Inputs Int -- int y - Global edge index. Note this is the unmarked y coordinate, -- i.e. the coordinate without scattering regions removed. -> Int -- int nx - Number of vertices left to right -> Int -- int ny - Number of vertices top to bottom -> Double -- float lx - Length of horizontal edges (distance between vertices in x direction) -> Double -- float ly - Length of vertical edges (distance between vertices in y direction) -> Double -- float k - Plane wave wavenumber -> Double -- float theta - Direction of wave propagation -> Double -- float E0 - Magnitude of incident plane wave -- Outputs -> Complex Double -- complex float e - Magnitude of electric field on edge y calcincidentfield y nx ny lx ly k theta e0 = let (xg, yg) = itoxy y nx ny in --Determine whether edge is horizontal or vertical let isvertical = xg >= nx in let (xvalueTmp, yvalueTmp) = edgetoxy y nx ny in let xvalue = xvalueTmp * lx in let yvalue = yvalueTmp * ly in -- Convert x and y edge coordinates to x and y values and caluculate field if isvertical then mkPolar (-cos(theta)*e0) ( -k*(xvalue*cos(theta)+yvalue*sin(theta))) else mkPolar (sin(theta)*e0) ( -k*(xvalue*cos(theta)+yvalue*sin(theta))) -- | Auxiliary function for oracle /b/. build_circuit getconnection :: Int -> Int -> Int -> Int -> Int -> Int getconnection y i' nx ny maxConnectivity = let i = getIntFromParam i' in let maxC = getIntFromParam maxConnectivity in let (ex,ey) = itoxy y nx ny in let x = if ( (ex < nx) && (ey /= 1) ) then case i' of 1 -> y-2*nx+1 2 -> y-nx 3 -> y-nx+1 4 -> y 5 -> y+nx-1 6 -> y+nx 7 -> y+2*nx-1 _ -> -1 else if ( (ex < nx) && (ey == 1) ) then case i' of 1 -> y 2 -> y+nx-1 3 -> y+nx 4 -> y+2*nx-1 _ -> -1 else if ( (ex >= nx) && (ex /= nx) && (ex /= 2*nx-1) ) then case i' of 1 -> y-nx 2 -> y-nx+1 3 -> y-1 4 -> y 5 -> y+1 6 -> y+nx-1 7 -> y+nx _ -> -1 else if ( (ex >= nx) && (ex == nx) ) then case i' of 1 -> y-nx+1 2 -> y 3 -> y+1 4 -> y+nx _ -> -1 else if ( (ex >= nx) && (ex == 2*nx-1) ) then case i' of 1 -> y-nx 2 -> y-1 3 -> y 4 -> y+nx-1 _ -> -1 else -1 in if (i > maxC) then -1 else if (x > nx*(ny-1)+ny*(nx-1)) then -1 else x -- | Auxiliary function to @'template_paramZero'@. build_circuit local_loop_with_index_aux :: Int -> Int -> t -> (Int -> t -> t) -> t local_loop_with_index_aux i n x f = case paramMinus n i of 0 -> x _ -> local_loop_with_index_aux (paramSucc i) n (f i x) f -- | Local version of @'loop_with_index'@, for lifting. build_circuit local_loop_with_index :: Int -> t -> (Int -> t -> t) -> t local_loop_with_index n x f = local_loop_with_index_aux paramZero n x f -- | Oracle /b/. build_circuit getKnownWeights :: Int -> Int -> Int -> [(Double,Double)] -> Double -> Double -> Double -> Double -> Double -> Int -> Complex Double getKnownWeights y nx ny scatteringnodes lx ly k theta e0 maxConnectivity = let makeConnections i connections = let x = getconnection y (paramSucc i) nx ny maxConnectivity in let t = not $ checkedge x scatteringnodes nx ny in (x,t):connections in let calcTang b (c,t) = if t then let matElt = calcmatrixelement y c nx ny lx ly k in let incField = calcincidentfield c nx ny lx ly k theta e0 in b - matElt * incField else b in let connections = local_loop_with_index maxConnectivity [] makeConnections in if (not $ checkedge y scatteringnodes nx ny) then 0.0 :+ 0.0 else foldl calcTang (0.0 :+ 0.0) connections ---------------------------------------------------------------------- ---------------------------------------------------------------------- -- Testing functions test_template_sinc = do f <- template_sinc r <- qinit (0 :: FDouble) f r return () test_template_itoxy = do f <- template_itoxy x <- qinit (0 :: FSignedInt) y <- qinit (0 :: FSignedInt) z <- qinit (0 :: FSignedInt) g <- f x h <- g y k <- h z return () test_template_edgetoxy = do f <- template_edgetoxy x <- qinit (0 :: FSignedInt) y <- qinit (0 :: FSignedInt) z <- qinit (0 :: FSignedInt) g <- f x h <- g y k <- h z return () test_template_calcRweights = do f <- template_calcRweights n1 <- qinit (0 :: FSignedInt) n2 <- qinit (0 :: FSignedInt) n3 <- qinit (0 :: FSignedInt) x1 <- qinit (0 :: FDouble) x2 <- qinit (0 :: FDouble) x3 <- qinit (0 :: FDouble) x4 <- qinit (0 :: FDouble) x5 <- qinit (0 :: FDouble) f1 <- f n1 f2 <- f1 n2 f3 <- f2 n3 g1 <- f3 x1 g2 <- g1 x2 g3 <- g2 x3 g4 <- g3 x4 g5 <- g4 x5 return () test_template_calcincidentfield = do y' <- qinit (0 :: FSignedInt) nx' <- qinit (0 :: FSignedInt) ny' <- qinit (0 :: FSignedInt) lx' <- qinit (0 :: FDouble) ly' <- qinit (0 :: FDouble) k' <- qinit (0 :: FDouble) theta' <- qinit (0 :: FDouble) e0' <- qinit (0 :: FDouble) f <- template_calcincidentfield f1 <- f y' f2 <- f1 nx' f3 <- f2 ny' f4 <- f3 lx' f5 <- f4 ly' f6 <- f5 k' f7 <- f6 theta' f8 <- f7 e0' return () test_template_calcmatrixelement = do y1 <- qinit (0 :: FSignedInt) y2 <- qinit (0 :: FSignedInt) nx <- qinit (0 :: FSignedInt) ny <- qinit (0 :: FSignedInt) lx <- qinit (0 :: FDouble) ly <- qinit (0 :: FDouble) k <- qinit (0 :: FDouble) f <- template_calcmatrixelement f1 <- f y1 f2 <- f1 y2 f3 <- f2 nx f4 <- f3 ny f5 <- f4 lx f6 <- f5 ly f7 <- f6 k return () test_template_getconnection = do y <- qinit (0 :: FSignedInt) nx <- qinit (0 :: FSignedInt) ny <- qinit (0 :: FSignedInt) f <- template_getconnection f1 <- f y f2 <- f1 6 f3 <- f2 nx f4 <- f3 ny f5 <- f4 7 return () test_template_checkedge = do y <- qinit (0 :: FSignedInt) nx <- qinit (0 :: FSignedInt) ny <- qinit (0 :: FSignedInt) s <- qinit [(0 :: FDouble,0 :: FDouble),(0 :: FDouble,0 :: FDouble)] f <- template_checkedge f1 <- f y f2 <- f1 s f3 <- f2 nx f4 <- f3 ny return ()