{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
{-# LINE 1 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
{-# 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).

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




{-# LINE 31 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| 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




 |] )
{-# LINE 32 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | Auxiliary function.

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)))



{-# LINE 40 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| 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)))



 |] )
{-# LINE 41 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | The function sin /x/ \/ /x/.

sinc :: Double -> Double
sinc x = if (x /= 0.0) then (sin x) / x else 1.0



{-# LINE 46 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| sinc :: Double -> Double
                      sinc x = if (x /= 0.0) then (sin x) / x else 1.0



 |] )
{-# LINE 47 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | Auxiliary function.

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)



{-# LINE 55 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| 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)



 |] )
{-# LINE 56 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | 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/).

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 lx/ly - k*k*lx*ly/3.0
                -- B12 and B21
                else if ( (abs(yg1-yg2) == 1) && (abs(xg1-xg2) == 0) && (xg1<nx) ) then -lx/ly + k*k*lx*ly/12.0
                -- B34 and B43
                else if ( (abs(yg1-yg2)==0) && (abs(xg1-xg2) == 1) && (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
                -- B41
                else if ( (yg2==(yg1+1)) && (xg2==(xg1-nx)) && (xg2<nx) ) then 1.0
                -- B42
                else if ( (yg1==yg2) && (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<nx) ) then 1.0
                -- B23
                else if ( (yg2==yg1) && (xg2==(xg1+nx-1)) && (xg1<nx) ) then 1.0
                else -1.0
        in
        let c = if ( (y1==y2) && ( (xg1==nx) || (xg1==(2*nx-1)) ) ) then 0.0 :+ (k*ly)
                else if  ( (y1==y2) && ( ((yg1==1) && (xg1<nx)) || ((yg1==ny) && (xg1<nx)) ) ) then 0.0 :+ (k*lx)
                else 0.0 :+ 0.0
        in (b :+ 0.0) + c



{-# LINE 118 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| calcmatrixelement ::

                              Int
                              -> Int
                              -> Int
                              -> Int
                              -> Double
                              -> Double
                              -> Double

                              -> Complex Double
                      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

                                      else if ( (y1==y2) && (xg1<nx) ) then lx/ly - k*k*lx*ly/3.0

                                      else if ( (abs(yg1-yg2) == 1) && (abs(xg1-xg2) == 0) && (xg1<nx) ) then -lx/ly + k*k*lx*ly/12.0

                                      else if ( (abs(yg1-yg2)==0) && (abs(xg1-xg2) == 1) && (xg1>=nx) ) then -ly/lx + k*k*lx*ly/12.0

                                      else if ( (yg1==(yg2+1)) && (xg1==(xg2-nx+1)) && (xg2>=nx) ) then -1.0

                                      else if ( (yg2==(yg1+1)) && (xg2==(xg1-nx+1)) && (xg1>=nx) ) then -1.0

                                      else if ( (yg1==(yg2+1)) && (xg1==(xg2-nx)) && (xg1<nx) ) then 1.0

                                      else if ( (yg2==(yg1+1)) && (xg2==(xg1-nx)) && (xg2<nx) ) then 1.0

                                      else if ( (yg1==yg2) && (xg1==(xg2+nx)) && (xg1>=nx) ) then -1.0

                                      else if ( (yg2==yg1) && (xg2==(xg1+nx)) && (xg2>=nx) ) then -1.0

                                      else if ( (yg1==yg2) && (xg1==(xg2+nx-1)) && (xg2<nx) ) then 1.0

                                      else if ( (yg2==yg1) && (xg2==(xg1+nx-1)) && (xg1<nx) ) then 1.0
                                      else -1.0
                              in
                              let c = if ( (y1==y2) && ( (xg1==nx) || (xg1==(2*nx-1)) ) ) then 0.0 :+ (k*ly)
                                      else if  ( (y1==y2) && ( ((yg1==1) && (xg1<nx)) || ((yg1==ny) && (xg1<nx)) ) ) then 0.0 :+ (k*lx)
                                      else 0.0 :+ 0.0
                              in (b :+ 0.0) + c



 |] )
{-# LINE 119 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | Auxiliary function.

get_edges l = case l of
                [] -> []
                (x:tt) -> case tt of
                   [] -> []
                   (y:t) -> (x,y):(get_edges (y:t))


{-# LINE 126 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| get_edges l = case l of
                                      [] -> []
                                      (x:tt) -> case tt of
                                         [] -> []
                                         (y:t) -> (x,y):(get_edges (y:t))


 |] )
{-# LINE 127 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | Auxiliary function.

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)





{-# LINE 138 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| 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)





 |] )
{-# LINE 139 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | Oracle /r/.

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<nx) ) then

         let i = (mkPolar lx (k*yc*sin(phi)))*
                 (mkPolar 1.0 (k*xc*cos(phi)))*
                 ((sinc (k*lx*(cos phi)/2.0)) :+ 0.0) in

         let r = ( (- sin phi) :+ k*ly )*((cos(theta - phi))/ly :+ 0.0) in i * r


     else if ( (yg==ny) && (xg<nx) ) then

         let i = (mkPolar lx (k*yc*sin(phi)))*
                 (mkPolar 1.0 (k*xc*cos(phi)))*
                 ((sinc (k*lx*(cos phi)/2.0)) :+ 0.0) in

         let r = ( (- sin phi) :+ (- k*ly) )*((cos(theta - phi)/ly) :+ 0.0) in i * r

     else 0.0 :+ 0.0




{-# LINE 185 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| 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<nx) ) then

                               let i = (mkPolar lx (k*yc*sin(phi)))*
                                       (mkPolar 1.0 (k*xc*cos(phi)))*
                                       ((sinc (k*lx*(cos phi)/2.0)) :+ 0.0) in

                               let r = ( (- sin phi) :+ k*ly )*((cos(theta - phi))/ly :+ 0.0) in i * r


                           else if ( (yg==ny) && (xg<nx) ) then

                               let i = (mkPolar lx (k*yc*sin(phi)))*
                                       (mkPolar 1.0 (k*xc*cos(phi)))*
                                       ((sinc (k*lx*(cos phi)/2.0)) :+ 0.0) in

                               let r = ( (- sin phi) :+ (- k*ly) )*((cos(theta - phi)/ly) :+ 0.0) in i * r

                           else 0.0 :+ 0.0




 |] )
{-# LINE 186 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | Auxiliary function for oracle /A/.

convertband :: Int -> 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





{-# LINE 238 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| convertband :: Int -> 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





 |] )
{-# LINE 239 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | 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/.

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)





{-# LINE 297 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| 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)





 |] )
{-# LINE 298 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | 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/.

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)))



{-# LINE 345 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| calcincidentfield ::

                                      Int

                              -> Int
                              -> Int
                              -> Double
                              -> Double
                              -> Double
                              -> Double
                              -> Double

                              -> Complex Double
                      calcincidentfield y nx ny lx ly k theta e0 =
                              let (xg, yg) = itoxy y nx ny in

                              let isvertical = xg >= nx in
                              let (xvalueTmp, yvalueTmp) = edgetoxy y nx ny in
                              let xvalue = xvalueTmp * lx in
                              let yvalue = yvalueTmp * ly in

                              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)))



 |] )
{-# LINE 346 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | Auxiliary function for oracle /b/.

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




{-# LINE 401 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| 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




 |] )
{-# LINE 402 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | Auxiliary function to @'template_paramZero'@.

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



{-# LINE 410 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| 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



 |] )
{-# LINE 411 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | Local version of @'loop_with_index'@, for lifting.

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



{-# LINE 416 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| 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



 |] )
{-# LINE 417 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
-- | Oracle /b/.

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






{-# LINE 438 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
$( decToCircMonad [d| 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






 |] )
{-# LINE 439 "Quipper/Algorithms/QLS/TemplateOracle.hs" #-}
----------------------------------------------------------------------
----------------------------------------------------------------------
-- 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 ()