text:
-- This file is part of Quipper. Copyright (C) 2011-2016. Please see the
-- file COPYRIGHT for a list of authors, copyright holders, licensing,
-- and other details. All rights reserved.
--
-- ======================================================================

{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeSynonymInstances #-}

{-# OPTIONS -fcontext-stack=50 #-}

-- | This module contains the Quipper implementation of the Quantum
-- Linear Systems Algorithm.
--
-- The algorithm estimates the radar cross section for a FEM
-- scattering problem by using amplitude estimation to calculate
-- probability amplitudes.
--
-- The notations are based on the paper
--
-- * B. D. Clader, B. C. Jacobs, C. R. Sprouse. Quantum algorithm to
-- calculate electromagnetic scattering cross
-- sections. <http://arxiv.org/abs/1301.2340>.
module Algorithms.QLS.QLS where

import Quipper

import QuipperLib.QFT
import QuipperLib.Arith
import QuipperLib.Decompose

import Data.Complex
import qualified Data.Map as Map

import qualified Algorithms.QLS.TemplateOracle as Template
import Algorithms.QLS.QDouble
import Algorithms.QLS.QSignedInt
import Algorithms.QLS.CircLiftingImport
import Algorithms.QLS.Utils

import Libraries.Auxiliary(boollist_of_int_bh)

-- | The type of 'oracle_A' input arguments during runtime.
type OracleARunTime = Double  -- ^Value resolution.
       -> Int                 -- ^Band.
       -> Bool                -- ^Argflag.
       -> ([Qubit],[Qubit],[Qubit])  -- ^(x=index,y+node,z+value).
       -> Circ ([Qubit],[Qubit],[Qubit])

-- | The type of 'oracle_b' and 'oracle_r' input arguments during runtime.
type OracleBRRunTime = Double  -- ^Magnitude resolution.
        -> Double              -- ^Phase resolution.
        -> ([Qubit],[Qubit],[Qubit]) -- ^(x=index,m+magnitude,p+phase).
        -> Circ ([Qubit],[Qubit],[Qubit])

-- | A type to encapsulate all three oracles.
data Oracle = Oracle {
  oracle_A :: RunTimeParam -> OracleARunTime,
  oracle_b :: RunTimeParam -> OracleBRRunTime,
  oracle_r :: RunTimeParam -> OracleBRRunTime
}

-- | A set of oracles using only blackboxes.
dummy_oracle :: Oracle
dummy_oracle = Oracle {
  oracle_A = \d r i b x -> named_gate "Oracle A" x,
  oracle_b = \d1 d2 r x -> named_gate "Oracle b" x,
  oracle_r = \d1 d2 r x -> named_gate "Oracle r" x
}

-- | A type to hold the runtime parameters.
data RunTimeParam = RT_param {
  k :: Double,       -- ^ Wave number.
  theta :: Double,   -- ^ Incident wave angle.
  phi :: Double,     -- ^ Direction of desired far field radiation pattern.
  e0 :: Double,      -- ^ Incident wave amplitude.
  lambda :: Double,  -- ^ Wavelength.
  xlength :: Double, -- ^ /x/-length of square scattering region.
  ylength :: Double, -- ^ /y/-length of square scattering region.

  scatteringnodes :: [(Int,Int)], -- ^ Metallic region.

  nx :: Int,
  ny :: Int,
  lx :: Double,
  ly :: Double,

  kappa :: Double,
  epsilon :: Double,
  t0 :: Double,
  r :: Double,
  b_max :: Double,
  r_max :: Double,
  d  :: Int,
  nb :: Int,
  p2 :: Double,
  n0 :: Int,
  n1 :: Int,
  n2 :: Int,
  n4 :: Int,

  -- argflags for oracle_A
  magnitudeArgflag :: Bool,
  phaseArgflag :: Bool
} deriving (Show)

-- | A convenient set of runtime parameters for testing.
dummy_RT_param :: RunTimeParam
dummy_RT_param = RT_param {

  k = 2.0*pi*1.0,
  theta = 0.0*pi/4.0,
  phi = 0.0,
  e0 = 1.0,
  lambda = (k dummy_RT_param)/(2.0*pi),
  xlength = 2.0*(lambda dummy_RT_param),
  ylength = 2.0*(lambda dummy_RT_param),

  scatteringnodes =
    let rt = dummy_RT_param in
    let xul = round((fromIntegral $ nx rt)/2)
              - round((xlength rt)/(2.0*(lx rt))) in -- Upper left x index
    let yul = round((fromIntegral $ ny rt)/2)
              - round((ylength rt)/(2.0*(ly rt))) in -- Upper left y index
    let xlr = round((fromIntegral $ nx rt)/2)
              + round((xlength rt)/(2.0*(lx rt))) in  -- Lower right x index
    let ylr = round((fromIntegral $ ny rt)/2)
              + round((ylength rt)/(2.0*(ly rt))) in -- Lower right y in
      [(xul, yul), (xlr, ylr)],

  nx = 12885,
  ny = 12885,
  lx = 0.1,
  ly = 0.1,

  kappa = 1.0,
  epsilon = 1.0,
  t0 = 1.0,
  r = 1.0,
  b_max = 1.0,
  r_max = 1.0,
  d  = 3,
  nb = 2,
  p2 = 3,
  n0 = 3,
  n1 = 3,
  n2 = 3,
  n4 = 3,

  -- argflags for oracle_A
  magnitudeArgflag = False,
  phaseArgflag = True
}

-- | A set of larger values, for testing scalability.
large_RT_param :: RunTimeParam
large_RT_param = RT_param {

  k = 2.0*pi*1.0,
  theta = 0.0*pi/4.0,
  phi = 0.0,
  e0 = 1.0,
  lambda = (k large_RT_param)/(2.0*pi),
  xlength = 2.0*(lambda large_RT_param),
  ylength = 2.0*(lambda large_RT_param),

  scatteringnodes =
    let rt = large_RT_param in
    let xul = round((fromIntegral $ nx rt)/2)
              - round((xlength rt)/(2.0*(lx rt))) in -- Upper left x index
    let yul = round((fromIntegral $ ny rt)/2)
              - round((ylength rt)/(2.0*(ly rt))) in -- Upper left y index
    let xlr = round((fromIntegral $ nx rt)/2)
              + round((xlength rt)/(2.0*(lx rt))) in  -- Lower right x index
    let ylr = round((fromIntegral $ ny rt)/2)
              + round((ylength rt)/(2.0*(ly rt))) in -- Lower right y in
      [(xul, yul), (xlr, ylr)],

  nx = 12885,
  ny = 12885,
  lx = 0.1,
  ly = 0.1,

  kappa = 1e4,
  epsilon = 0.01,
  t0 = 7.0e6,
  r = 2.5e12,
  b_max = 5.0,
  r_max = 1.01,
  d  = 7,
  nb = 9,
  p2 = (  1.0 / (4-(4**(1/3)))  ),
  n0 = 14,
  n1 = 24,
  n2 = 30,
  n4 = 65,

  -- argflags for oracle_A
  magnitudeArgflag = False,
  phaseArgflag = True
}

-- | A set of smaller values, for manageable yet meaningful output.
small_RT_param :: RunTimeParam
small_RT_param = RT_param {

  k = 2.0*pi*1.0,
  theta = 0.0*pi/4.0,
  phi = 0.0,
  e0 = 1.0,
  lambda = (k small_RT_param)/(2.0*pi),
  xlength = 2.0*(lambda small_RT_param),
  ylength = 2.0*(lambda small_RT_param),

  scatteringnodes =
    let xul = 2 in
    let yul = 2 in
    let xlr = 3 in
    let ylr = 3 in
      [(xul, yul), (xlr, ylr)],

  nx = 4,
  ny = 4,
  lx = 0.1,
  ly = 0.1,

  kappa = 1e4,
  epsilon = 0.01,
  t0 = 7.0e6,
  r = 2.5e12,
  b_max = 5.0,
  r_max = 1.01,
  d  = 7,
  nb = 9,
  p2 = (  1.0 / (4-(4**(1/3)))  ),
  n0 = 14,
  n1 = 24,
  n2 = 6,
  n4 = 65,

  -- argflags for oracle_A
  magnitudeArgflag = False,
  phaseArgflag = True
}

-- | Apply an [exp −/iYt/] gate. The timestep /t/ is a parameter.
expYt :: Timestep -> Qubit -> Circ Qubit
expYt = named_rotation "exp(-i%Y)"

-- | Apply an [exp −/iYt/] gate. The timestep /t/ is a parameter.
expYt_at :: Timestep -> Qubit -> Circ ()
expYt_at = named_rotation_at "exp(-i%Y)"

-- | Read a list of bits and make it into a 'Double', by multiplying
-- its integer value by the provided factor.
dynamic_lift_double :: Double -> [Bit] -> Circ Double
dynamic_lift_double factor cl = do
      cdiscard cl

-- Implementation note: removed dynamic_lift as for now it breaks all
-- output formats except ASCII.

--      bl <- dynamic_lift cl
      let sign = 1 -- if (head $ reverse bl) then 1 else -1
      let unsigned_value = 1 -- integer_of_intm_unsigned $
                              -- intm_of_boollist_bh (tail $ reverse bl)
      return (sign * factor * (fromIntegral unsigned_value))

-- | A black box gate to stand in as a replacement for QFT.
qft_for_show :: [Qubit] -> Circ [Qubit]
qft_for_show qs = named_gate "QFT" qs

-- | Main function: for estimating the radar cross section for a FEM
-- scattering problem. The problem can be reduced to the calculation
-- of four angles: φ[sub /b/], φ[sub /bx/], φ[sub /r/1] and φ[sub /r/0].
qlsa_FEM_main :: RunTimeParam -> Oracle -> Circ Double
qlsa_FEM_main param oracle = do
     comment "FEM_main"
     phi_b  <- qlsa_AmpEst_phi_b param oracle
     phi_bx <- qlsa_AmpEst_phi_bx param oracle
     phi_r1 <- qlsa_AmpEst_phi_bxr param oracle True
     phi_r0 <- qlsa_AmpEst_phi_bxr param oracle False
     let sigma = ((((fromIntegral $ nb param) ^ 2) * ((b_max param) ^ 2) * ((r_max param) ^ 2) * ((sin phi_b) ^ 2)) / ( 4 * pi))
     comment "FEM_main"
     return sigma

-- * Amplitude Estimation Functions

-- | Estimates φ[sub /b/], related to the probability of success for the
-- preparation of the known state /b/, using amplitude amplification.
qlsa_AmpEst_phi_b :: RunTimeParam -> Oracle -> Circ Double
qlsa_AmpEst_phi_b param oracle = do
    g <- qinit $ replicate (n0 param) False
    with_ancilla_init (replicate (n2 param) False) $ \x -> do
       label (g,x) ("g","x")
       with_ancilla $ \a -> do
           label (a) ("anc. a")
           with_ancilla $ \b -> do
               label (b) ("anc. b")
               g <- map_hadamard g
               u_b (x,b)
               loop g u_g (x,b,a)
               return ()
           return ()
    g' <- qft_big_endian g -- QFT : Is it really big-endian?
    value_bits  <- measure g'
    value_double <- dynamic_lift_double 1.0 value_bits
    return (pi * value_double / (2 ** (fromIntegral $ n0 param)))
    where
        loop :: [Qubit] -> (a -> Circ ()) -> a -> Circ ()
        loop [] f x = return ()
        loop (h:t) f x = do
            f x `controlled` h
            loop t f' x;
            where
               f' x = do f x; f x

        u_b :: ([Qubit],Qubit) -> Circ ()
        u_b xb = qlsa_StatePrep param xb (oracle_b oracle param) (1.0/(b_max param))

        u_g :: ([Qubit],Qubit,Qubit) -> Circ ()
        u_g (x,b,a) = do
            comment "U_g starts"
            gate_Z_at b
            -- For unitrary linear transformation, adjoint == inverse, hence reverse_...
            (reverse_generic_imp u_b) (x,b)
            qnot_at a `controlled` x .==. (map (\x -> False) x)
            gate_X_at a
            gate_Z_at a
            gate_X_at a
            qnot_at a `controlled` x .==. (map (\x -> False) x)
            u_b (x,b)
            comment "U_g ends"
            return ()

-- | Testing function for 'qlsa_AmpEst_phi_b'.
test_qlsa_AmpEst_phi_b :: Bool -> IO ()
test_qlsa_AmpEst_phi_b dummyRTParamFlag = do
    let param = if dummyRTParamFlag then dummy_RT_param else large_RT_param
    print_simple GateCount (qlsa_AmpEst_phi_b param dummy_oracle)
    print_simple Preview (qlsa_AmpEst_phi_b param dummy_oracle)

-- | Estimates φ[sub /bx/], related to the probability of success in
-- computing solution value /x/.
qlsa_AmpEst_phi_bx :: RunTimeParam -> Oracle -> Circ Double
qlsa_AmpEst_phi_bx param oracle  = do
    g <- qinit (take (n0 param) (repeat False))
    with_ancilla_init (take (n2 param) (repeat False)) $ \x -> do
      with_ancilla $ \a -> do
          with_ancilla $ \b -> do
              with_ancilla $ \s -> do
                  g <- map_hadamard g
                  u_bx (x,b,s)
                  loop g u_g (x,b,s,a)
                  return ()
              return ()
          return ()
    g' <- qft_big_endian g -- QFT : Is it really big-endian?
    value_bits  <- measure g'
    value_double <- dynamic_lift_double 1.0 value_bits
    return (pi * value_double / (2 ** (fromIntegral $ n0 param)))
    where
        loop :: [Qubit] -> (a -> Circ ()) -> a -> Circ ()
        loop [] f x = return ()
        loop (h:t) f x = do
            f x `controlled` h
            loop t f' x
            where
            f' x = do f x; f x

        u_bx :: ([Qubit],Qubit,Qubit) -> Circ ()
        u_bx (x,b,s) = do
            qlsa_StatePrep param (x,b) (oracle_b oracle param) (1.0/(b_max param))
            qlsa_Solve_x param (x,s) oracle
            return ()

        u_g :: ([Qubit],Qubit,Qubit,Qubit) -> Circ ()
        u_g (x,b,s,a) = do --named_gate_at "Ug_phi_bx" (x,b,s,a)
            qnot_at a `controlled` b .&&. s
            gate_Z_at a
            qnot_at a `controlled` b .&&. s
            (reverse_generic_imp u_bx) (x,b,s)
            qnot_at a `controlled` [ q .==. 0 | q <- x ] .&&. b .==. 0 .&&. s .==. 0
            gate_X_at a
            gate_Z_at a
            gate_X_at a
            qnot_at a `controlled` [ q .==. 0 | q <- x ] .&&. b .==. 0 .&&. s .==. 0
            u_bx (x,b,s)
            return ()

-- | Estimates φ[sub /r/0] and φ[sub /r/1] (depending on the boolean
-- parameter), related to the overlap of the solution with the
-- arbitrary state /r/.
qlsa_AmpEst_phi_bxr :: RunTimeParam -> Oracle -> Bool -> Circ Double
qlsa_AmpEst_phi_bxr param oracle target = do
    g <- qinit (take (n0 param) (repeat False))
    with_ancilla_init (take (n2 param) (repeat False)) $ \x -> do
      with_ancilla_init (take (n2 param) (repeat False)) $ \y -> do
        with_ancilla $ \a -> do
          with_ancilla $ \b -> do
             with_ancilla $ \s -> do
                with_ancilla $ \r -> do
                    with_ancilla $ \c -> do
                        g <- map_hadamard g
                        u_r (x,y,b,s,r,c)
                        loop g u_g (x,y,b,s,r,c,a)
                        return ()
                    return ()
                return ()
             return ()
          return ()
    g' <- qft_big_endian g -- QFT : Is it really big-endian?
    value_bits  <- measure g'
    value_double <- dynamic_lift_double 1.0 value_bits
    return (pi * value_double / (2 ** (fromIntegral $ n0 param)))
    where
    loop :: [Qubit] -> (a -> Circ ()) -> a -> Circ ()
    loop [] f x = return ()
    loop (h:t) f x = do
       f x `controlled` h
       loop t f' x
       where
         f' x = do f x; f x

    u_r :: ([Qubit],[Qubit],Qubit,Qubit,Qubit,Qubit) -> Circ ()
    u_r (x,y,b,s,r,c) = do
        qlsa_Solve_xr param (x,y,b,s,r,c) oracle
        return ()

    u_g :: ([Qubit],[Qubit],Qubit,Qubit,Qubit,Qubit,Qubit) -> Circ ()
    u_g (x,y,b,s,r,c,a) = do --named_gate_at "Ug_phi_bxr" (x,y,b,s,r,c,a)
         qnot_at a `controlled` (b .==. 1 .&&. s .==. 1 .&&. r .==. 1 .&&. c .==. target)
         gate_Z_at a
         qnot_at a `controlled` (b .==. 1 .&&. s .==. 1 .&&. r .==. 1 .&&. c .==. target)
         (reverse_generic_imp u_r) (x,y,b,s,r,c)
         qnot_at a `controlled` [ q .==. 0 | q <- x ] .&&. [ q .==. 0 | q <- y ] .&&. b .==. 0 .&&. s .==. 0 .&&. r .==. 0 .&&. c .==. 0
         gate_X_at a
         gate_Z_at a
         gate_X_at a
         qnot_at a `controlled` [ q .==. 0 | q <- x ] .&&. [ q .==. 0 | q <- y ] .&&. b .==. 0 .&&. s .==. 0 .&&. r .==. 0 .&&. c .==. 0
         u_r (x,y,b,s,r,c)
         return ()

-- * State Preparation.

-- | Prepares a quantum state /x/, as specified by an oracle function,
-- entangled with a single qubit flag /q/ marking the desired state.
qlsa_StatePrep ::
    RunTimeParam
    -> ([Qubit], Qubit) -- x & qare handles to wires to be changed by StatePrep
    -> OracleBRRunTime -- Common type of (oracle_b oracle) and (oracle_r oracle), if oracle is typed Oracle
    -> Double
    -> Circ ()
qlsa_StatePrep param (x, q) oracle phi0 = do
  _ <- (flip $ box ("qlsa_StatePrep_" ++ (show phi0))) (x,q) $ \(x,q) -> do
    -- comment "StatePrep starts"
    label (x,q) ("x", "q")
    with_ancilla_list (n4 param) $ \m -> do
        label (m) ("anc. m")
        with_ancilla_list (n4 param) $ \p -> do
            label (p) ("anc. p")
            x <- map_hadamard x
            (x, m, p) <- oracle phi0 (epsilon param) (x, m, p)
            qlsa_ControlledPhase p (epsilon param) False
            qlsa_ControlledRotation (m, q) phi0 False
            (x, m, p) <- oracle phi0 (epsilon param) (x, m, p)
            return (x,q)
    -- comment "StatePrep ends"
  return ()

-- | Testing function for 'qlsa_StatePrep'.
test_qlsa_StatePrep :: Bool -> IO ()
test_qlsa_StatePrep dummyRTParamFlag = do
          let param = if dummyRTParamFlag then dummy_RT_param else large_RT_param
          let oraclebRunTime = (oracle_b dummy_oracle param)
          let testCirc = do
              x <- qinit $ replicate (n2 param) False
              q <- qinit False
              qlsa_StatePrep param (x,q) oraclebRunTime 1.0
          print_simple GateCount testCirc
          print_simple Preview testCirc

-- * Linear System Solver Functions

-- | Implements the QLSA procedure to generate the solution state |/x/〉.
qlsa_Solve_x :: RunTimeParam ->([Qubit],Qubit) -> Oracle -> Circ ()
qlsa_Solve_x param (x,s) oracle = do
  _ <- (flip $ box "qlsa_Solve_x") (x,s) $ \(x,s) -> do
    with_ancilla_list (n1 param) $ \t -> do
        with_ancilla_list (length t) $ \f -> do
            let phi0 = 2 * pi  / (2 ** (fromIntegral $ n1 param) - (epsilon param))
            t <- map_hadamard t
            u_hs (t,x)
            t <- qft_big_endian t
            integer_inverse (t,f)
            qlsa_ControlledRotation (f,s) phi0 False
            integer_inverse (t,f)
            (reverse_generic_endo qft_big_endian) t
            (reverse_generic_imp u_hs) (t,x)
            t <- map_hadamard t
            return()
        return ()
    return (x,s)
  return ()
    where
        u_hs :: ([Qubit],[Qubit]) -> Circ ()
        u_hs (t,x) = do
            qlsa_HamiltonianSimulation param (t,x) (oracle_A oracle param)
            return ()

-- | Implementation of the integer division. The two registers are
-- supposed to be of the same size and represent little-headian
-- unsigned integers, i.e., the head of the list holds the least
-- significant bit.
integer_inverse :: ([Qubit],[Qubit]) -> Circ ()
integer_inverse (t,f) = do
  _ <- (flip $ box "integer_inverse") (t,f) $ \(t,f) -> do
    -- sanity check
    if (length t /= length f)
      then error "integer_inverse: registers of distinct sizes"
      else return ()
    -- initialize an unsigned integer to 2^(length t) - 1
    with_ancilla_init (map (\_ -> True) t) $ \num -> do
      -- perform the division (encapsulated in a subroutine)
      let d = classical_to_reversible $ \(t,num) -> do
                let x = ((qdint_of_qulist_lh num),(qdint_of_qulist_lh t))
                (_,_,f') <- uncurry q_div_unsigned x
                return $ qulist_of_qdint_lh f'
      d ((t,num),f)
      return (t,f)
  return ()

-- | Implements the complete QLSA procedure to find the
-- solution state |/x/〉 and then implements the swap protocol
-- required for estimation of 〈/x/|/r/〉.
qlsa_Solve_xr :: RunTimeParam ->([Qubit],[Qubit],Qubit,Qubit,Qubit,Qubit) -> Oracle -> Circ ()
-- qlsa_Solve_xr param (x,y,b,s,r,c) oracle  = named_gate_at "Solve_xr" (x,y,b,s,r,c)
qlsa_Solve_xr param (x,y,b,s,r,c) oracle = do
  _ <- (flip $ box "qlsa_Solve_xr") (x,y,b,s,r,c) $ \(x,y,b,s,r,c) -> do
    qlsa_StatePrep param (x,b) (oracle_b oracle param) (1.0/(b_max param))
    qlsa_Solve_x param (x,s) oracle
    qlsa_StatePrep param (y,r) (oracle_r oracle param) (1.0/(r_max param))
    hadamard_at c
    swap_at y x  `controlled` c
    hadamard_at c
    return (x,y,b,s,r,c)
  return ()

-- * Hamiltonian Simulation Functions.

-- | Uses a quantum register |/t/〉 to control the
-- implementation of the Suzuki method for simulating a Hamiltonian
-- specified by an oracle function.
qlsa_HamiltonianSimulation ::
    RunTimeParam
    -> ([Qubit], [Qubit])
    -> OracleARunTime
    -> Circ ()
qlsa_HamiltonianSimulation param (t, x) oracleA = do
  _ <- (flip $ box "qlsa_HamiltonianSimulation") (t,x) $ \(t,x) -> do
    -- Code the first line in a way that depends on the length of t rather
    -- explicitly on (n1 param) which is supposed to be the length of t
    -- Hence replaced the line below by the line after it.
    -- let denom = 2 * (r param) * (fromIntegral ((n1 param) - 1))
    let denom = 2 * (r param) * ( 2^((length t) - 1) )
    let t1 = (p2 param) * (t0 param) / denom
    let t2 = (1.0 - 4.0 * (p2 param)) * (t0 param) / denom
    (t,x) <- box_loopM  "TrotterLoop" (round $ r param) (t,x) (hs_loop t1 t2)
    return (t,x)
  return ()
    where
          hs_loop :: Double -> Double -> ([Qubit], [Qubit]) -> Circ ([Qubit], [Qubit])
          hs_loop t1 t2 (t,x) = do
            u_z_at (t,x) t1
            u_z_at (t,x) t1
            u_z_at (t,x) t2
            u_z_at (t,x) t1
            u_z_at (t,x) t1
            return (t,x)

          u_z_at :: ([Qubit], [Qubit]) -> Double -> Circ ()
          u_z_at (t, x) timeStep = do
              for (nb param) 1 (-1) $ \jj -> do
                  qlsa_HsimKernel param (t, x) jj timeStep oracleA
                  endfor
              for 1 (nb param) 1 $ \jj -> do
                  qlsa_HsimKernel param (t, x) jj timeStep oracleA
                  endfor
              return ()

-- | Testing function for 'qlsa_HamiltonianSimulation'.
test_qlsa_HamiltonianSimulation :: Bool -> IO ()
test_qlsa_HamiltonianSimulation dummyRTParamFlag = do
    let param = if dummyRTParamFlag then dummy_RT_param else large_RT_param
    let oracleARunTime = (oracle_A dummy_oracle param)
    let testCirc = do
              t <- qinit (take (n0 param) (repeat False))
              x <- qinit (take (n4 param) (repeat False))
              label(t,x) ("t", "x")
              qlsa_HamiltonianSimulation param (t,x) oracleARunTime
    print_simple GateCount testCirc
--    print_simple Preview testCirc

-- | Uses an oracle function and timestep control register (/t/) to
-- apply 1-sparse Hamiltonian to the input state |/t/, /x/〉.
qlsa_HsimKernel ::
    RunTimeParam
    -> ([Qubit], [Qubit])
    -> Int
    -> Double
    -> OracleARunTime
    -> Circ ()
-- qlsa_HsimKernel param tx band timeStep oracleA = named_gate_at "HsimKernel" tx

qlsa_HsimKernel param (t, x) band timeStep oracleA = do
  _ <- (flip $ box ("qlsa_HsimKernel" ++ (show band) ++ (show timeStep))) (t,x) $ \(t,x) -> do
    let phiP = (epsilon param)
    with_ancilla_list (n2 param) $ \y -> do
       with_ancilla_list (n4 param) $ \m -> do
           with_ancilla_list (n4 param) $ \p -> do
               label (y,m,p) ("y","m","p")
               -- phases
               oracleA phiP band (phaseArgflag param) (x,y,p)
               qlsa_ControlledPhase p phiP False
               oracleA phiP band (phaseArgflag param) (x,y,p)
               -- magnitudes
               let phiMag = 2 ** (negate $ fromIntegral $ after_radix_length)
               oracleA phiMag band (magnitudeArgflag param) (x,y,m)
               for 0 ((length t) - 1) 1 $ \ii -> do
                   let phi_mt = timeStep * phiMag * (2^ii)
                   qlsa_ApplyHmag param (x,y,m) phi_mt `controlled` (t !! ii)
                   endfor
               oracleA phiMag band (magnitudeArgflag param) (x,y,m)
               -- phases again
               oracleA phiP band (phaseArgflag param) (x,y,p)
               qlsa_ControlledPhase p phiP True
               oracleA phiP band (phaseArgflag param) (x,y,p)
               return ()
           return ()
       return ()
    return (t,x)
  return ()

-- | Testing function for 'qlsa_HsimKernel'.
test_qlsa_HsimKernel :: Bool -> IO ()
test_qlsa_HsimKernel dummyRTParamFlag = do
    let param = if dummyRTParamFlag then dummy_RT_param else large_RT_param
    let oracleARunTime = (oracle_A dummy_oracle param)
    let testCirc = do
              t <- qinit (take (n0 param) (repeat False))
              x <- qinit (take (n4 param) (repeat False))
              label(t,x) ("t", "x")
              qlsa_HsimKernel param (t,x) 1 0.1 oracleARunTime
    print_simple GateCount testCirc
--    print_simple Preview testCirc

-- | Applies the magnitude component of coupling elements in a
-- 1-sparse Hamiltonian.
qlsa_ApplyHmag :: RunTimeParam -> ([Qubit], [Qubit], [Qubit]) -> Double -> Circ ()
qlsa_ApplyHmag param (x,y,m) phi0 = do
  _ <- (flip $ box ("qlsa_ApplyHmag " ++ (show phi0))) (x,y,m) $ \(x,y,m) -> do
    let (onOne, onZero) = (True, False)
    with_ancilla $ \a -> do -- Assume initialized to False (0)
        label (a) ("anc. a")
        if (length x /= length y)
           then error "qlsa_ApplyHmag: Input registers x and y have different lengths."
           else return ()
        let length_xy = length x
        for 0 (length_xy - 1) 1 $ \ii -> do
            let (xi, yi) = (x !! ii, y !! ii)
            w (xi, yi)
            qnot_at a `controlled` (xi .==. onOne .&&. yi .==. onZero)
            endfor
        qlsa_ControlledPhase (m ++ [a]) phi0 False
        for (length_xy - 1) 0 (-1) $ \ii -> do
            let (xi, yi) = (x !! ii, y !! ii)
            qnot_at a `controlled` (xi .==. onOne .&&. yi .==. onZero)
            w (xi, yi)
            endfor
        return ()
    return (x,y,m)
  return ()

-- | Testing function for 'qlsa_ApplyHmag'.
test_qlsa_ApplyHmag :: Bool -> IO ()
test_qlsa_ApplyHmag dummyRTParamFlag = do
    let param = if dummyRTParamFlag then dummy_RT_param else large_RT_param
    let testCirc = do
              x <- qinit (take (n2 param) (repeat False))
              y <- qinit (take (n2 param) (repeat False))
              m <- qinit (take (n4 param) (repeat False))
              label(x,y,m) ("x", "y", "m")
              qlsa_ApplyHmag param (x,y,m) 0.1
    print_simple GateCount testCirc
    print_simple Preview testCirc

-- | Auxiliary function: the /W/-gate.
w :: (Qubit, Qubit) -> Circ ()
-- w xy = named_gate_at "W" xy
w (xi,yi) = do
  _ <- box "w" (\(xi, yi) -> do
    label (xi,yi) ("x[i]","y[i]")
    gate_X_at xi `controlled` yi
    gate_X_at yi `controlled` xi
    hadamard_at yi `controlled` xi
    gate_X_at yi `controlled` xi
    gate_X_at xi `controlled` yi
    return (xi,yi)) (xi,yi)
  return ()

-- | Testing function for 'w'.
test_w :: IO ()
test_w = do
    let testCirc = do
        xi <- qinit False
        yi <- qinit False
        w (xi, yi)
    print_simple GateCount testCirc
    print_simple Preview testCirc

-- * Controlled Logic Operations

-- | Applies a phase shift of φ\/2 to the signed input register |φ〉.

-- For c, bit locations are counted starting from 0 (least significant) to (length c - 2) (most significant)
-- last c (or c !! (n-1)) is the sign bit
qlsa_ControlledPhase :: [Qubit] -> Double -> Bool -> Circ ()
-- qlsa_ControlledPhase c phi0 f = named_gate_at "CPhase" c
qlsa_ControlledPhase c phi0 f = do
  _ <- (flip $ box ("qlsa_ControlledPhase " ++ (show phi0) ++ " " ++ (show f))) c $ \c -> do
    with_ancilla $ \a -> do -- ancilla is initialized to False
        if f then (qnot_at a) else return ()
        let signQubit = last c
        qnot_at a `controlled` signQubit
        for 0 (length c - 2) 1 $ \ii -> do
            let theta = ( (2.0^(ii)) * phi0 / 2.0) -- Note the divide by 2
            -- The following three lines are equivalent
            expZt_at theta a `controlled` (c !! ii)
            -- a <- expZt theta a `controlled` (c !! ii); return ()
            -- qlsa_ControlledU (c !! ii, a) (expZt theta) False
            endfor
        qnot_at a `controlled` signQubit
        return c
  return ()

-- | Applies a rotation of φ\/2 to the signed input register |φ〉.
qlsa_ControlledRotation :: ([Qubit], Qubit) -> Double -> Bool -> Circ ()
-- qlsa_ControlledRotation ct phi0 f = named_gate_at "CRotate" ct
qlsa_ControlledRotation (c, t) phi0 f = do
  _ <- (flip $ box ("qlsa_ControlledRotation " ++ (show phi0) ++ " " ++ (show f))) (c,t) $ \(c,t) -> do
    if f then (qnot_at t) else return ()
    let signQubit = last c
    qnot_at t `controlled` signQubit
    for 0 (length c - 2) 1 $ \ii -> do
        let theta = (2.0^(ii)) * phi0 / 2.0
        expYt_at theta t `controlled` (c !! ii)
        endfor
    qnot_at t `controlled` signQubit
    if f then (qnot_at t) else return ()
    return (c,t)
  return ()

----------------------------------------------------------------------
-- * Oracles

-- | Map a 'QDouble' into an integer, understood as being scaled by
-- the given factor. Take the factor and the size of the output
-- register as parameter.
make_factor_rep :: Double -> Int -> QDouble -> Circ [Qubit]
make_factor_rep factor size p = do
     -- number of high bits of p
     let p_int_size = (xdouble_length p) - (xdouble_exponent p)
     -- get the required number of bits for doing the encoding
     let auxsize = max (size - 1) p_int_size
     -- do the operation:
     --    (1) build a QDouble from the factor: need (size-1) bits of integer part
     qfactor <- qinit $ fdouble_of_double (xdouble_exponent p)
                                          (auxsize + xdouble_exponent p) factor
     --    (2) scale p
     p_large <- qdouble_pad_to_extent (auxsize,- (xdouble_exponent p)) p
     --    (3) divide p by factor
     qreal_multiples <- unpack template_symb_slash_ p_large qfactor
     --    (4) get the floor of the result
     q_floor <- template_floor
     qmultiples <- q_floor qreal_multiples
     --    (5) set them in the right format (it has size elements)
     let (SInt tp bp) = qmultiples
     let new_p = (take (size - 1) $ reverse tp) ++ [bp]
     return new_p

-- | Implements the oracle for the arbitrary state |/r/〉, using the
-- Template Haskell implementation of 'calcRweights'.
inline_oracle_r :: RunTimeParam -> Double -> Double -> ([Qubit],[Qubit],[Qubit]) -> Circ ([Qubit],[Qubit],[Qubit])
inline_oracle_r rt factor_m factor_p =  box ("Or " ++ (show factor_m) ++ " " ++ (show factor_p)) $ decompose_generic Toffoli $ \(x',m',p') ->
    with_ancilla_init (fromIntegral $ nx rt :: FSignedInt) $ \qnx ->
     with_ancilla_init (fromIntegral $ ny rt :: FSignedInt) $ \qny ->
      with_ancilla_init (fdouble $ lx rt) $ \qlx ->
       with_ancilla_init (fdouble $ ly rt) $ \qly ->
        with_ancilla_init (fdouble $ k rt) $ \qk ->
         with_ancilla_init (fdouble $ theta rt) $ \qtheta ->
          with_ancilla_init (fdouble $ phi rt) $ \qphi ->
                -- the x register is smaller than the size of a QSInt,
                -- and x is an unsigned integer: make up a positive
                -- sign and a pad

                with_ancilla_init False $ \bx -> do
                  with_ancilla_init (replicate (fixed_int_register_length - (n2 rt))
                                               False) $ \pad_x -> do

                    let x = SInt (pad_x ++ (reverse x')) bx
                    let f = classical_to_reversible $
                              \(x,nx,ny,lx,ly,k,theta,phi) -> do
                                 (m,p) <- unpack Template.template_calcRweights
                                            x nx ny lx ly k theta phi
                                 new_p <- make_factor_rep factor_p (n4 rt) p
                                 new_m <- make_factor_rep factor_m (n4 rt) m
                                 return (new_m, new_p)
                    f ((x,qnx,qny,qlx,qly,qk,qtheta,qphi),(m',p'))
                    return (x',m',p')

-- | Implements the oracle for the known state |/b/〉, using the
-- Template Haskell implementation of 'getKnownWeights'.
inline_oracle_b :: RunTimeParam -> Double -> Double -> ([Qubit],[Qubit],[Qubit]) -> Circ ([Qubit],[Qubit],[Qubit])
inline_oracle_b rt factor_m factor_p = box ("Ob " ++ (show factor_m) ++ " " ++ (show factor_p)) $ decompose_generic Toffoli $ \(x',m',p') ->
    -- Make ancillas for constant values

    with_ancilla_init (listpair_fmap fromIntegral $ scatteringnodes rt :: [(FDouble,FDouble)]) $ \qscatteringnodes ->
     with_ancilla_init (fromIntegral $ ny rt :: FSignedInt) $ \qny ->
      with_ancilla_init (fdouble $ lx rt) $ \qlx ->
       with_ancilla_init (fdouble $ ly rt) $ \qly ->
        with_ancilla_init (fdouble $ k rt) $ \qk ->
         with_ancilla_init (fdouble $ theta rt) $ \qtheta ->
          with_ancilla_init (fdouble $ e0 rt) $ \qe0 ->
           with_ancilla_init (fromIntegral $ nx rt :: FSignedInt) $ \qnx -> do

                -- the x register is smaller than the size of a QSInt,
                -- and x is an unsigned integer: make up a positive
                -- sign and a pad

                with_ancilla_init False $ \bx -> do
                  with_ancilla_init (replicate (fixed_int_register_length - (n2 rt))
                                               False) $ \pad_x -> do

                    let x = SInt (pad_x ++ (reverse x')) bx

                    let f = classical_to_reversible $
                              \(y,nx,ny,scatteringnodes,lx,ly,k,theta,e0) -> do
                                  -- get some QSInt and QDouble
                                  (m,p) <- unpack Template.template_getKnownWeights
                                              y nx ny scatteringnodes lx ly k theta e0 7
                                  new_p <- make_factor_rep factor_p (n4 rt) p
                                  new_m <- make_factor_rep factor_m (n4 rt) m

                                  return (new_m, new_p)

                    f ((x,qnx,qny,qscatteringnodes,qlx,qly,qk,qtheta,qe0),(m',p'))
                    return (x',m',p')

-- | Implementation of the oracle calculating the matrix /A/
-- corresponding to the discretization of the scattering problem,
-- using the Template Haskell implementation of
-- 'getNodeValuesMoreOutputs'.
inline_oracle_A ::  RunTimeParam -> Double -> Int -> Bool -> ([Qubit],[Qubit],[Qubit]) -> Circ ([Qubit],[Qubit],[Qubit])
inline_oracle_A rt factor band argflag =  box ("OA " ++ (show band) ++ " " ++ (show argflag)) $ decompose_generic Toffoli $ \(x',y',p') -> do

    let argflag' = if argflag then PTrue else PFalse

    -- Make ancillas for constant values

    with_ancilla_init (listpair_fmap fromIntegral $ scatteringnodes rt :: [(FDouble, FDouble)]) $ \qscatteringnodes ->
      with_ancilla_init (fromIntegral $ ny rt :: FSignedInt) $ \qny ->
        with_ancilla_init (fromRational $ toRational $ lx rt) $ \qlx ->
          with_ancilla_init (fromRational $ toRational $ ly rt) $ \qly ->
            with_ancilla_init (fromRational $ toRational $ k rt) $ \qk ->
             with_ancilla_init (fromIntegral $ nx rt :: FSignedInt) $ \qnx -> do

                -- the x register is smaller than the size of a QSInt,
                -- and x is an unsigned integer: make up a positive
                -- sign and a pad

                with_ancilla_init False $ \bx -> do
                  with_ancilla_init (replicate (fixed_int_register_length - (n2 rt))
                                               False) $ \pad_x -> do

                    let x = SInt (pad_x ++ (reverse x')) bx

                    let f = classical_to_reversible $
                              \(v,nx,ny,scatteringnodes,lx,ly,k) -> do
                               -- get some QSInt and QDouble
                               (y,p) <- unpack Template.template_getNodeValuesMoreOutputs
                                           v band nx ny scatteringnodes lx ly k argflag' 7

                               new_p <- make_factor_rep factor (n4 rt) p

                               -- set y in the correct format (it has n2 elements)
                               -- we can assume that the sign is positive.
                               -- we also assume that n2 < size of QSInt register
                               let (SInt ty _) = y
                               let new_y = (take (n2 rt) $ reverse ty)

                               -- return the values in the right format.
                               return (new_y,new_p)

                    f ((x,qnx,qny,qscatteringnodes,qlx,qly,qk),(y',p'))

                    return (x',y',p')

-- | Encapsulate the inline oracles in Template Haskell into an object
-- of type 'Oracle'.
inline_oracle :: Oracle
inline_oracle = Oracle {
  oracle_A = inline_oracle_A,
  oracle_b = inline_oracle_b,
  oracle_r = inline_oracle_r
}

max_stars_repo_name: silky/quipper

max_issues_repo_name: qbit-/quipper

max_forks_repo_name: qbit-/quipper


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# OPTIONS_GHC -fno-warn-unused-imports #-}

-- | A custom prelude
--
-- $license

module StdLib (

      identity
    , (&)
    , uncons
    , applyN
    , print
    , io

    , LText
    , LByteString

    , CIText
    , CIByteString

      -- * List

    , (L.\\)
    , head
    , nubOrd
    , nubOn
    , sortOn
    , count
    , chunksOf

      -- * Show

    , Print(..)
    , putText
    , putLText
    , tshow
    , bshow
    , blshow

      -- * Bool

    , bool
    , whenM
    , unlessM
    , ifM
    , guardM
    , while
    , whileM

      -- * Debug

    , undefined
    , error
    , trace
    , traceM
    , traceIO
    , traceShow
    , traceShowM
    , notImplemented

      -- * Monad

    , Monad(..)
    , MonadPlus(..)

    , (=<<)
    , (>=>)
    , (<=<)
    , forever

    , join
    , mfilter
    , filterM
    , mapAndUnzipM
    , zipWithM
    , zipWithM_
    , foldM
    , foldM_
    , replicateM
    , replicateM_
    , concatMapM
    , concatForM
    , whenJust
    , unlessNull

    , guard
    , when
    , unless

    , liftM
    , liftM2
    , liftM3
    , liftM4
    , liftM5
    , liftM'
    , liftM2'
    , ap

    , (<$!>)

      -- * Applicative

    , orAlt
    , orEmpty
    , eitherA

      -- * Either

    , maybeToLeft
    , maybeToRight
    , leftToMaybe
    , rightToMaybe
    , eitherToMaybe

      -- * Text

    , txtpack
    , txtunpack
    , ltxtpack
    , ltxtunpack
    , txt2ltxt
    , txt2bs
    , bs2txt
    , lbs2txt
    , ltxt2lbs
    , bs2ltxt
    , lbs2bs
    , bs2lbs

      -- * Set

    , setFromList

      -- * Map

    , mapFromList

      -- * Case-insensitive text

    , ci
    , ciOriginal

      -- * Time

    , Milliseconds
    , utcTimeToEpochTime

      -- * Async

    , async_

      -- * Re-exports

    , (P.$)
    , (P.$!)
    , P.minBound
    , P.maxBound

#if __GLASGOW_HASKELL__ >= 808
    , MonadFail(..)
#endif

    , module X
    ) where

import qualified Control.Monad.IO.Class as IO (liftIO)
import qualified Data.ByteString as X (ByteString)
import qualified Data.ByteString as BS
import qualified Data.ByteString.Char8 as BS8
import qualified Data.ByteString.Lazy
import qualified Data.ByteString.Lazy as BL
import qualified Data.ByteString.Lazy.Char8 as BL8
import qualified Data.CaseInsensitive as CI
import qualified Data.List as L
import qualified Data.Map as Map
import qualified Data.Set as Set
import qualified Data.Text as T
import qualified Data.Text.Encoding as T
import qualified Data.Text.IO as T
import qualified Data.Text.Lazy
import qualified Data.Text.Lazy as TL
import qualified Data.Text.Lazy.Encoding as TL
import qualified Data.Text.Lazy.IO as TL
import qualified Debug.Trace as T
import qualified Prelude as P
import qualified Prelude

import           Control.Applicative as X (Alternative (..), Applicative (..),
                                           Const (..), ZipList (..), liftA,
                                           liftA2, liftA3, optional, (<**>))
import           Control.Concurrent as X hiding (throwTo)
import           Control.Concurrent.Async as Async
import           Control.Monad hiding ((<$!>))
import           Control.Monad.Except as X (Except, ExceptT, MonadError,
                                            catchError, runExcept, runExceptT,
                                            throwError)
import           Control.Monad.IO.Class as X (MonadIO)
import           Control.Monad.Reader as X (MonadReader, Reader, ReaderT, ask,
                                            asks, local, runReader, runReaderT)
import           Control.Monad.ST as X
import           Control.Monad.State as X (MonadState, State, StateT,
                                           evalStateT, execState, execStateT,
                                           get, gets, modify, put, runStateT,
                                           withState)
import           Control.Monad.Trans as X (MonadTrans (..), lift)
import           Control.Monad.Trans.Maybe as X (runMaybeT)
import           Data.Bits as X
import           Data.Bool as X hiding (bool)
import           Data.ByteString as X (ByteString)
import           Data.Char as X (Char, isSpace)
import           Data.Complex as X
import           Data.Either as X
import           Data.Eq as X
import           Data.Foldable as X hiding (foldl1, foldr1, maximum, maximumBy,
                                     minimum, minimumBy)
import           Data.Function as X (const, fix, flip, on, (.))
import           Data.Functor as X (Functor (..), void, ($>), (<$>))
import           Data.Functor.Identity as X
import           Data.Int as X
import           Data.IntMap as X (IntMap)
import           Data.IntSet as X (IntSet)
import           Data.List as X (break, drop, filter, intercalate, isPrefixOf,
                                 replicate, reverse, sortBy, splitAt, take)
import           Data.List.NonEmpty as X (NonEmpty (..), nonEmpty)
import           Data.Map as X (Map)
import           Data.Maybe as X hiding (fromJust)
import           Data.Monoid as X hiding ((<>))
import           Data.Ord as X
import           Data.Semigroup as X hiding (First, Last, getFirst, getLast)
import           Data.Sequence as X (Seq)
import           Data.Set as X (Set)
import           Data.Text as X (Text)
import           Data.Time as X
import           Data.Time.Clock.POSIX
import           Data.Traversable as X
import           Data.Tuple as X
import           Data.Word as X
import           GHC.Enum as X (Enum (..))
import           GHC.Exts as X (Constraint, FunPtr, Ptr, the)
import           GHC.Float as X hiding (log)
import           GHC.Int as X
import           GHC.IO as X (IO)
import           GHC.Num as X
import           GHC.Real as X
import           GHC.Show as X
import           GHC.Stack
import           Numeric as X (showHex)
import           System.Directory as X (createDirectoryIfMissing, doesFileExist,
                                        getCurrentDirectory,
                                        setCurrentDirectory)
import           System.Environment as X (getArgs, getEnvironment, lookupEnv,
                                          setEnv, unsetEnv)
import           System.Exit as X
import           System.FilePath as X (FilePath, takeBaseName, (</>))
import           System.IO as X (BufferMode (..), Handle, hClose, hFlush,
                                 hSetBinaryMode, hWaitForInput)
import           System.Process as X (readProcess, readProcessWithExitCode)
import           Text.Printf as X (PrintfArg, hPrintf, printf)
import           Text.Read as X (Read, readEither, readMaybe, reads)

type LText = Data.Text.Lazy.Text

type LByteString = Data.ByteString.Lazy.ByteString

type CIText = CI.CI Text

type CIByteString = CI.CI ByteString

infixl 1 &

(&) :: a -> (a -> b) -> b
x & f = f x

identity :: a -> a
identity x = x

uncons :: [a] -> Maybe (a, [a])
uncons []     = Nothing
uncons (x:xs) = Just (x, xs)

applyN :: Int -> (a -> a) -> a -> a
applyN n f = X.foldr (.) identity (X.replicate n f)

print :: (MonadIO m, Show a) => a -> m ()
print = io . putStrLn . tshow

io :: MonadIO m => IO a -> m a
io = IO.liftIO

-- * List

head :: Foldable f => f a -> Maybe a
head = foldr (\x _ -> return x) Nothing

sortOn :: Ord o => (a -> o) -> [a] -> [a]
sortOn = sortBy . comparing

-- O(n * log n)
nubOrd :: Ord a => [a] -> [a]
nubOrd l = go Set.empty l
  where
    go _ []     = []
    go s (x:xs) =
      if x `Set.member` s
      then go s xs
      else x : go (Set.insert x s) xs

nubOn :: Eq b => (a -> b) -> [a] -> [a]
nubOn f = fmap snd
        . L.nubBy ((==) `on` fst)
        . fmap (\x -> let y = f x in y `Prelude.seq` (y, x))

count :: Eq a => a -> [a] -> Int
count x = length . filter (x ==)

chunksOf :: Int -> [a] -> [[a]]
chunksOf _ [] = []
chunksOf n xs = take n xs : chunksOf n (drop n xs)

-- * Show

class Print a where
    putStr :: MonadIO m => a -> m ()
    putStrLn :: MonadIO m => a -> m ()

instance Print Text where
    putStr = io . T.putStr
    putStrLn = io . T.putStrLn

instance Print LText where
    putStr = io . TL.putStr
    putStrLn = io . TL.putStrLn

instance Print BS.ByteString where
    putStr = io . BS.putStr
    putStrLn = io . BS8.putStrLn

instance Print LByteString where
    putStr = io . BL.putStr
    putStrLn = io . BL8.putStrLn

instance Print [Char] where
    putStr = io . Prelude.putStr
    putStrLn = io . Prelude.putStrLn

putText :: MonadIO m => Text -> m ()
putText = putStrLn
{-# SPECIALIZE putText :: Text -> IO () #-}

putLText :: MonadIO m => LText -> m ()
putLText = putStrLn
{-# SPECIALIZE putLText :: LText -> IO () #-}

tshow :: Show a => a -> Text
tshow = T.pack . show

{-# WARNING bshow "'bshow' remains in code" #-}
bshow :: Show a => a -> ByteString
bshow = BS8.pack . show

{-# WARNING blshow "'blshow' remains in code" #-}
blshow :: Show a => a -> LByteString
blshow = BL8.pack . show

-- * Bool

bool :: a -> a -> Bool -> a
bool f t p = if p then t else f

whenM :: Monad m => m Bool -> m () -> m ()
whenM p m = p >>= flip when m

unlessM :: Monad m => m Bool -> m () -> m ()
unlessM p m = p >>= flip unless m

ifM :: Monad m => m Bool -> m a -> m a -> m a
ifM p x y = p >>= \b -> if b then x else y

guardM :: MonadPlus m => m Bool -> m ()
guardM f = guard =<< f

while :: Monad m => m Bool -> m a -> m ()
while cond action = do
    c <- cond
    when c P.$ do
        _ <- action
        while cond action

-- | While loop. run 'f' while 'p'
whileM :: (MonadPlus m1, Monad m) => (a -> m Bool) -> m a -> m (m1 a)
whileM p f = go
  where
    go = do
      x <- f
      r <- p x
      if r then (return x `mplus`) `liftM` go else return mzero

{-# WARNING undefined "'undefined' remains in code" #-}
undefined :: HasCallStack => a
undefined = P.undefined

{-# WARNING error "'error' remains in code" #-}
error :: HasCallStack => P.String -> a
error = P.error

{-# WARNING trace "'trace' remains in code" #-}
trace :: P.String -> a -> a
trace = T.trace

{-# WARNING traceShow "'traceShow' remains in code" #-}
traceShow :: Show a => a -> a
traceShow a = T.trace (P.show a) a

{-# WARNING traceShowM "'traceShowM' remains in code" #-}
traceShowM :: (Show a, Monad m) => a -> m ()
traceShowM a = T.traceM (P.show a)

{-# WARNING traceM "'traceM' remains in code" #-}
traceM :: Monad m => P.String -> m ()
traceM = T.traceM

{-# WARNING traceIO "'traceIO' remains in code" #-}
traceIO :: P.String -> IO ()
traceIO = T.traceIO

{-# WARNING notImplemented "'notImplemented' remains in code" #-}
notImplemented :: a
notImplemented = P.error "Not implemented"

-- * Monad

concatMapM :: Monad m => (a -> m [b]) -> [a] -> m [b]
concatMapM f xs = liftM concat (mapM f xs)

whenJust :: Monad m => Maybe a -> m () -> m ()
whenJust = when . isJust

unlessNull :: Monad m => [a] -> m () -> m ()
unlessNull = unless . L.null

concatForM :: (Monad m, Traversable t) => t a -> (a -> m [b]) -> m [b]
concatForM ls m = fmap concat (forM ls m)

liftM' :: Monad m => (a -> b) -> m a -> m b
liftM' = (<$!>)
{-# INLINE liftM' #-}

liftM2' :: Monad m => (a -> b -> c) -> m a -> m b -> m c
liftM2' f a b = do
    x <- a
    y <- b
    let z = f x y
    z `Prelude.seq` return z
{-# INLINE liftM2' #-}

(<$!>) :: Monad m => (a -> b) -> m a -> m b
f <$!> m = do
    x <- m
    let z = f x
    z `Prelude.seq` return z
{-# INLINE (<$!>) #-}

-- * Applicative

orAlt :: (Alternative f, Monoid a) => f a -> f a
orAlt f = f <|> pure mempty

orEmpty :: Alternative f => Bool -> a -> f a
orEmpty b a = if b then pure a else empty

eitherA :: Alternative f => f a -> f b -> f (Either a b)
eitherA a b = (Left <$> a) <|> (Right <$> b)

-- * Either

leftToMaybe :: Either l r -> Maybe l
leftToMaybe = either Just (const Nothing)

rightToMaybe :: Either l r -> Maybe r
rightToMaybe = either (const Nothing) Just

maybeToRight :: l -> Maybe r -> Either l r
maybeToRight l = maybe (Left l) Right

maybeToLeft :: r -> Maybe l -> Either l r
maybeToLeft r = maybe (Right r) Left

eitherToMaybe :: Either a b -> Maybe b
eitherToMaybe (Left _)  = Nothing
eitherToMaybe (Right r) = Just r

-- * Text / ByteString

txtpack :: Prelude.String -> Text
txtpack = T.pack

txtunpack :: Text -> Prelude.String
txtunpack = T.unpack

ltxtpack :: Prelude.String -> LText
ltxtpack = TL.pack

ltxtunpack :: LText -> Prelude.String
ltxtunpack = TL.unpack

txt2ltxt :: Text -> LText
txt2ltxt = TL.fromStrict

txt2bs :: Text -> ByteString
txt2bs = T.encodeUtf8

bs2txt :: ByteString -> Text
bs2txt = T.decodeUtf8

lbs2txt :: LByteString -> Text
lbs2txt = TL.toStrict . TL.decodeUtf8

ltxt2lbs :: LText -> LByteString
ltxt2lbs = TL.encodeUtf8

bs2ltxt :: LByteString -> LText
bs2ltxt = TL.decodeUtf8

lbs2bs :: LByteString -> ByteString
lbs2bs = BL.toStrict

bs2lbs :: ByteString -> LByteString
bs2lbs = BL.fromChunks . (:[])

-- * Set

setFromList :: Ord a => [a] -> Set a
setFromList = Set.fromList

-- * Map

mapFromList :: Ord a => [(a, b)] -> Map a b
mapFromList = Map.fromList

-- * Case-insensitive

ci :: CI.FoldCase s => s -> CI.CI s
ci = CI.mk

ciOriginal :: CI.CI a -> a
ciOriginal = CI.original

-- * Time

type Milliseconds = Double

utcTimeToEpochTime :: UTCTime -> Integer
utcTimeToEpochTime = round . utcTimeToPOSIXSeconds

-- * Async

-- | Just like 'Async.async' but ignores the return value 'Async'
--
-- An 'Async' is a Haskell RTS thread that writes to a TMVar when
-- finished.
--
-- 'cancel' will send the Haskell thread a @SIGKILL@ signal but you don't
-- need to explicitly kill the thread. If the 'Async' can be garbage
-- collected, the thread will still run to completion and then be
-- properly cleaned up.
--
-- However, if there is an exception in the 'Async', then an exception
-- won't be propagate to the "parent" thread (unless 'wait' is somehow
-- used)
async_ :: IO a -> IO ()
async_ = void . Async.async

ext: hs


text:
{-# LANGUAGE MultiParamTypeClasses #-}
module ChaosBox.Geometry.Class
  ( Boundary(..)
  , HasP2(..)
  , Intersects(..)
  , getP2
  , setP2
  , modifyP2
  )
where

import           ChaosBox.Geometry.P2
import           Control.Lens         (Lens', lens, (%~), (&), (.~), (^.))
import           Data.Complex
import           Linear.V2

-- | Class of objects that can be queried for points
class Boundary a where
  containsPoint :: a -> P2 -> Bool

class HasP2 a where
  _V2 :: Lens' a P2

instance HasP2 P2 where
  _V2 = _xy

instance HasP2 (Complex Double) where
  _V2 = lens (\(a :+ b) -> V2 a b) (\_ (V2 x y) -> x :+ y)

getP2 :: HasP2 a => a -> P2
getP2 = (^. _V2)

setP2 :: HasP2 a => a -> P2 -> a
setP2 x p2 = x & _V2 .~ p2

modifyP2 :: HasP2 a => a -> (P2 -> P2) -> a
modifyP2 x f = x & _V2 %~ f

class Intersects a b where
  intersectionPoints :: a -> b -> [P2]
  intersects :: a -> b -> Bool
  intersects a b = null (intersectionPoints a b)

max_stars_repo_name: akrmn/chaosbox

max_issues_repo_name: akrmn/chaosbox

max_forks_repo_name: akrmn/chaosbox


text:
{-# LANGUAGE ViewPatterns #-}

module Main where

import Codec.Image.DevIL
import Data.Array.Unboxed
import Math.Probably.MCMC
import Math.Probably.Sampler
import Math.Probably.RandIO
import Math.Probably.FoldingStats
import Control.Monad.State.Strict
import System.Cmd
import System.Environment
import Data.Array.IO
import System.IO

import Data.Array.Unboxed
import Data.Array.Base
import qualified Data.Vector.Unboxed as V
import qualified Data.StorableVector as SV
import qualified Data.StorableVector.Base as SVB
import Numeric.LinearAlgebra hiding (find)
import qualified Math.Probably.PDF as PDF
import Data.List
import Data.Maybe
import Control.Applicative
import Data.Ord
import System.Directory

import Data.IORef
import System.IO.Unsafe

import CVUtils
import Edge
import CVhelp
import Foreign.C.String
import Foreign.C.Types
import Unsafe.Coerce

bgImRef :: IORef Image
bgImRef = unsafePerformIO $ newIORef undefined

visibleMat :: IORef BitImage
visibleMat = unsafePerformIO $ newIORef undefined

nogoMat :: IORef BitImage
nogoMat = unsafePerformIO $ newIORef undefined

inNogoMat x y = unsafePerformIO $ do
    nm <- readIORef nogoMat
    return $ readBitImage nm x y

triRef :: IORef Triangles
triRef = unsafePerformIO $ newIORef undefined

gloTris = unsafePerformIO $ readIORef triRef

sdv = 0.5
sdrot = 0.4
sddiv = 2
sddivRot = 4
sdlen = 0.1
sdSideDisp = 0.3

hiddenflat = 1

nparticles = 1000

gauss2i :: R -> R-> R-> R-> Sampler (R,R)
gauss2i m1 sd1 m2 sd2 = do
   u1 <- unitSample
   u2 <- unitSample
   let l = sqrt(-2*log(u1))
       th = 2*pi*u2
--   return (sqrt(-2*log(u1))*cos(2*pi*u2)*sd1+m1,
--           sqrt(-2*log(u1))*sin(2*pi*u2)*sd2+m2)
   return (l*cos(th)*sd1+m1,
           l*sin(th)*sd2+m2)

evolve :: Int -> Obj -> Sampler (Obj,Obj)
evolve retries o@(Obj vlen side len x y rot) = do
        (nvlen, nside) <- gauss2i vlen sdv 0 sdSideDisp
--        nside <- gaussD 0 sdSideDisp
        (nrot, nlen) <- gauss2i rot sdrot len sdlen
--        nlen <- gaussD len sdlen
        let no = Obj nvlen nside nlen
                                 (x+nvlen*cos nrot+nside*cos(nrot-pi/2))
                                 (y+nvlen*sin nrot+nside*sin(nrot-pi/2)) nrot
        if inNogoMat (round $ posx no) (round $ posy no) && retries > 0
           then evolve (retries-1) o
           else return $ (o,no)

invisible :: [((Int,Int),(Int,Int))]

invisible = [((362,621),(355,654)),
             ((442,649),(443,674)),
             ((866,626),(870,658)),
             ((995,628),(999,657)),
             ((1149,629),(1167,658)),
             ((1003,71),(1010,36)),
             ((422,71),(417,42)),
             ((157,162),(127,158)),
             ((1178,469),(1187,501))]

--[ ((963,626),
--               (999,663))]

visible :: Int -> Int -> Bool
visible x y  = not $ any (\sqr -> within sqr  x  y) invisible

mkVisibleMat :: IO ()
mkVisibleMat = do
    tris <- loadTriangles
    writeIORef triRef tris
    let visIm = mkBitImage 1280 720
                 $ \x -> \y-> any (pointInTriangle (fromList [realToFrac x,
                                                              realToFrac y]))
                                  tris
                              && visible x y
    writeIORef visibleMat visIm
    let nogoIm = mkBitImage 1280 720
                 $ \x -> \y-> any (\sqr-> within sqr (realToFrac x) (realToFrac y)) nogo
                              || (not $ any (pointInTriangle (fromList [realToFrac x, realToFrac y])) tris)
    writeIORef nogoMat nogoIm
    return ()

objCentreInvisible :: Obj -> Bool
objCentreInvisible (Obj _ _ _ cx cy _) = not $ visible (round cx) (round cy)

nogo::  [Square]
nogo =
    [((1128,657.5),(1162,627.5)),
     ((963.25,656),(997,626.5)),
     ((836.5,653.75),(866.75,624.5)),
     ((389.75,619.5),(359.75,648.75)),
     ((162.75,309),(132,278.5)),
     ((158.5,191.25),(127.25,160.75)),
     ((422,72.25),(450.75,40.5)),
     ((757.5,72),(728.5,41.75)),
     ((1002.75,70.5),(974.25,39.75)),
     ((1191,194),(1221,222.75)),
     ((1185,435.25),(1213.5,463.5))]
--  [ ((963,626),
--     (997,663))]

nogoPrior :: Obj -> R
nogoPrior o = sum $ map f nogo
          where  f sqr | within sqr (posx o) (posy o) = -1e100
                       | otherwise = 0

{-nogoObj :: Obj -> Bool
nogoObj o = any (\sqr-> within sqr (posx o) (posy o)) nogo
              || (not $ any (pointInTriangle (fromList [posx o,
                                                        posy o])) gloTris) -}

prevprior :: Obj -> (Obj -> R)
prevprior (Obj vlen side len x y rot)
      (Obj nvlen nside nlen _ _ nrot)
        =   PDF.gaussD vlen (sdv/sddiv) nvlen -- "heavy tailed proposals"
          + PDF.gaussD rot (sdrot/sddivRot) nrot
          + PDF.gaussD len (sdlen/sddiv) nlen
          + PDF.gaussD 0 (sdSideDisp/sddiv) nside
          + PDF.gaussD 5 0.5 nlen
          + uniformLogPdf 3 7 nlen

track :: StaticParams -> String -> Int -> Int -> Int -> Obj -> StateT Seed IO Obj
track sp vidfnm startfrm nframes frameOffset obj0 = do
  let outFnm = fileroot vidfnm ++ ".pos"
  h<- lift $ openFile outFnm WriteMode
  lift $ withCString vidfnm open_video
  res <-  go h (replicate nparticles obj0) [startfrm..startfrm+nframes-1] obj0
  lift $ hClose h
  lift $ close_video
  return res
   where
    go _ objs [] o = return o
    go h objs (i:is) _ = do
           bgIm <- lift $ readIORef bgImRef
           nogomat <- lift $ readIORef nogoMat
--           lift $ putStrLn $"~/cvutils/extract "++vidfnm++" "++show i
--           lift $ system $"~/cvutils/extract "++vidfnm++" "++show i
--           frame <- lift $ readImage "extract.png"

           lift advance
           diffObjs <- sample $ mapM (evolve 10) objs
           --lift $ print2  "pixel1Red = " (frame!(0,0,0))
           vismat <- lift $ readIORef visibleMat
           let wparticles = {-# SCC "wpart" #-} (dropLosers $ particleLike sp bgIm vismat nogomat diffObjs)
           let anyInvisible = {-# SCC "anyvis" #-}any objCentreInvisible $ map fst wparticles
           let npart = if anyInvisible then 5*nparticles else nparticles
           let wparticles' = if anyInvisible then map (\(x,w) -> (x,w/hiddenflat)) wparticles else wparticles

           let smws = sumWeights wparticles'
           let cummSmws = cummWeightedSamples wparticles'
           lift $ print2 "fr,npart=" (i,npart)
           --lift $ print2 "winners= " (map snd wparticles')
           --lift $ hFlush stdout
           nextObjs <- {-# SCC "resample" #-} sample $ sequence
                              $ replicate npart
                              $ do u <- unitSample
                                   return . fst . fromJust $ find ((>=u*smws) . snd) cummSmws
           --lift $ print $ map snd $ take 10 $ reverse $ cummSmws

           let mobj = {-# SCC "mobj" #-}
                   runStat (pure Obj <*> before meanF vellen
                                     <*> before meanF sideDisp
                                     <*> before meanF objlen
                                     <*> before meanF posx
                                     <*> before meanF posy
                                     <*> before meanF rot)
                              nextObjs
           --lift $ putStrLn $ show (i,mobj, snd $ last $  wparticles)
           lift $ hPutStrLn h $ show (i,mobj, snd $ last $  wparticles)
           lift $ hFlush h
           when (i `rem` 50 == 0) $ {-# SCC "upbg" #-} do
             lift $ updateBgIm mobj
             when (i `rem` 100 == 0) $ {-# SCC "outframe" #-} do
               lift $ system $"~/cvutils/extract "++vidfnm++" "++show i
               frame <- lift $ readImage "extract.png"
               bgIm2 <- lift $ readIORef bgImRef
               lift $ writeImage ("bgtrack"++show (i+frameOffset)++".png") bgIm2
               markedIm1 <- lift $ markEllipse sp (mobj) frame
               markedIm <- lift $ markObjsOnImage nextObjs markedIm1
               lift $ writeImage ("frame"++show (i+frameOffset)++".png") markedIm
           go h nextObjs is mobj
           --return () -- $ mobj:rest

cu2w8 :: CUChar -> Word8
cu2w8 = unsafeCoerce

pixval :: Int -> Int -> CUChar -> Word8
pixval x y c = cu2w8 $ pixel_value (fromIntegral x) (fromIntegral (719-y)) c

particleLike :: StaticParams -> Image -> BitImage -> BitImage ->[(Obj,Obj)] -> [(Obj,R)]
particleLike sp@(SP noise len ecc) bgim vismat nogomat objprs = map pL objprs where
  radiusi = 2* ceiling (len*ecc) + 2
  objs = map snd objprs
  xmin = minOn (posx) objs - radiusi
  xmax = maxOn (posx) objs+radiusi
  ymin = minOn (posy) objs-radiusi
  ymax = maxOn posy objs+radiusi
  xys = [(x,y) | x <- [xmin.. xmax], y <- [ymin..ymax]]
  yrng = (ymax-ymin+1)
--  cvec = fromList [20,25,30]
  ifvis x y n = if readBitImage vismat x y then n else 0
  ifInside, ifoutside :: UArray (Int,Int) R
  ifInside = {-# SCC "ifinside" #-} listArray ((xmin,ymin),(xmax,ymax))
               [ifvis x y $ (gaussRnn noise 0.07 $ pixval x y 2)
                           + (gaussRnn noise 0.09 $ pixval x y 1)
                           + (gaussRnn noise 0.12 $ pixval x y 0)
                       | (x,y) <- xys]
  ifoutside = {-# SCC "ifoutside" #-} listArray ((xmin,ymin),(xmax,ymax))
              [ifvis x y $ (gaussW8nn noise (bgim!(y,x,0)) $ pixval x y 2)
                           + (gaussW8nn noise (bgim!(y,x,1)) $ pixval x y 1)
                           + (gaussW8nn noise (bgim!(y,x,2)) $ pixval x y 0)
                       | (x,y) <- xys]

  pL (old,o@(Obj  _ _ len cx cy rot)) = {-# SCC "pL" #-}
    let lenecc=len*ecc
        srot = lenecc*sin rot
        crot = lenecc*cos rot
        f1x = cx+crot
        f1y = cy+srot
        f2x = cx-crot
        f2y = cy-srot
--        f :: (R,R) -> R
        f sm xy@(x, y)
--          | not $ readBitImage vismat x y  = 0
          = {-# SCC "f" #-} if dist  f1x f1y   x  y  + dist  f2x f2y   x  y  < 2 * len
                               then sm+ ifInside!xy
                               else sm+ifoutside!xy
--                               then sm + (V.unsafeIndex) ifInside  ((x-xmin)*yrng+(y-ymin))
 --                              else sm + (V.unsafeIndex) ifoutside ((x-xmin)*yrng+(y-ymin))
--                               then (SVB.unsafeIndex) ifInside
--                               else (SVB.unsafeIndex) ifoutside ((x-xmin)*yrng+(y-ymin))

    in if readBitImage nogomat (round $ cx) (round $ cy)
          then (o,-1e100)
          else (o, prevprior old o + {-# SCC "fsum" #-} (noise * foldl' f 0 xys))
                   --x <- xs,
                   --y <- ys]))

updateBgIm (Obj _ _ _ cx cy _) = do
    im <- readIORef bgImRef
    mutIm <- thaw im
    ((loy,lox,_),(hiy,hix,_)) <- getBounds mutIm
    forM_ [(y,x,c) | x<- [lox..hix],
                     y<- [loy..hiy],
                     c<- [0..2],
                     dist cx cy x y > 30] $ \ix@(y,x,c)-> do
      now <- readArray mutIm ix
      writeArray mutIm ix $ (now `div` 2)  + (pixval x y (rgb2bgr c) `div` 2)
    newbg <- freeze (mutIm::MImage)
    writeIORef bgImRef newbg

rgb2bgr :: Int -> CUChar
rgb2bgr 2 = 0
rgb2bgr 0 = 2
rgb2bgr 1 = 1

{-testCV fvid = do
     system $ "~/cvutils/extract "++fvid++" 0"
     frame <- readImage "extract.png"
     print2 "via png " $ frame!(100,100,0)
     withCString fvid open_video
     advance
     print2 "via c " $ pixel_value 100 (719-100) 2
     print2 "via c " $ pixval 100 100 2
     close_video -}
     --error "foo"

main = do
     ilInit
     bgnm : fvid : (read -> x) : (read -> y) : (read -> rot) : rest <- getArgs
     --testCV fvid

     bgIm <-readImage bgnm
     writeIORef bgImRef bgIm
--     marked <- markBg bgIm
--     writeImage "markbg.png" marked
     mkVisibleMat
     print fvid
     vm <- readIORef visibleMat
     writeVisMat "vismat.png" vm bgIm
     let (frStart,nfr) = case rest of
            [] -> (0,6000)
            s:[] -> (read s, 6000)
            s1:s2:_ -> (read s1,read s2)

     --frame0 <-readImage "extract.png"
     {-marked <- markObjsOnImage [Obj (x,y) 0,
                                Obj (x+1,y) 0,
                                Obj (x,y+1) 0,
                                Obj (x+1,y+1) 0] frame0
     writeImage "marked.png" marked -}
     let --posterior = posteriorV bgIm frame0 (round x,round y)
         --postAndV v = (posteriorV bgIm frame0 (round x,round y) $ fromList v, v)
         sp = (SP 218 5 0.9)
--         rot = negate (pi/7)
         --initialsV = fromList [x,y,218, 6, rot, 0.9]
         initObj =  (Obj 1 rot 6 x y rot)
     --ellim <- markEllipse sp initObj frame0
     --writeImage "markell.png" ellim
     --print $ nogoPrior (Obj 0 rot 992 650 rot)
     --print $ nogoPrior (Obj 0 rot 997 650 rot)
     ----print $ visible 994 650
     --print $ visible 1000 650

     runRIO $ do
--         iniampar <- sample $ initialAdaMet 500 1e-4 posterior initialsV
--         AMPar v _ _ _ <- runAndDiscard 5000 (show . ampPar) iniampar $ adaMet False posterior
--         lift $ print v
--         track (v@> 2) (v @> 3) bgIm fvid 3 (v@>0,v@>1)
         if "%d" `isInfixOf` fvid
            then trackMany sp fvid initObj
            else track sp fvid frStart nfr 0 $ initObj

     return ()

putS s  = lift . putStrLn $ s

trackMany sp fvid initObj = go 0 0 initObj where
  go n frOff obj = do
     let fnm = takeWhile (/='%') fvid ++ show n ++ ".avi"
     putS $ "looking for "++fnm
     exists <- lift $ doesFileExist fnm
     if exists then reallyGo n fnm frOff obj else putS "all done" >> return obj
  reallyGo n fnm frOff obj0= do
     lift $ system $ "~/cvutils/extract "++fnm++" "++show i ++">tempout"
     lns <- lift $ lines `fmap` readFile "tempout"
     let nframes = read $ tail $ dropWhile (/='=') $ head $ filter ("total frames=" `isPrefixOf`) lns
     objNext <- track sp fnm 0 (nframes-2) frOff $ obj0
     go (n+1) (frOff+nframes-2) objNext

-- initial on wl0: (956,641)

--rm marked.png && track mixed.png ~/Dropbox/woodlice/wl0.avi '(956,641)' && eog marked.png
--rm -f frame*.png && sudo cabal install --global && track mixed.png ~/Dropbox/woodlice/wl0.avi 956 641 '-0.448'

--rm -f frame*.png && sudo cabal install --global && track mixed.png ~/Dropbox/woodlice/wl0.avi 1092 518 '-5.12'

max_stars_repo_name: glutamate/cvutils

max_issues_repo_name: glutamate/cvutils

max_forks_repo_name: glutamate/cvutils


text:
{- |
Copyright: (c) 2020 Thomas Tuegel
SPDX-License-Identifier: BSD-3-Clause
Maintainer: Thomas Tuegel <ttuegel@mailbox.org>

-}

module Injection
    ( Injection (..)
    , Retraction (..)
    ) where

import Data.Complex (Complex ((:+)))
import Data.Dynamic (Dynamic, Typeable, fromDynamic, toDyn)
import Data.Fixed (Fixed, HasResolution)
import Data.Functor.Const (Const (..))
import Data.Functor.Identity (Identity (..))
import Data.List.NonEmpty (NonEmpty (..))
import Data.Maybe (maybeToList)
import Data.Monoid (Dual (..))
import Data.Monoid (Product (..))
import Data.Monoid (Sum (..))
import Data.Monoid (Any (..))
import Data.Monoid (All (..))
import qualified Data.Monoid as Monoid (First (..), Last (..))
import Data.Ord (Down (..))
import Data.Ratio (Ratio)
import qualified Data.Ratio as Ratio
import Data.Semigroup (Max (..), Min (..))
import qualified Data.Semigroup as Semigroup (First (..), Last (..))
import Data.Text (Text)
import qualified Data.Text as Text
import qualified Data.Text.Lazy as Lazy (Text)
import qualified Data.Text.Lazy as Text.Lazy
import Data.Void (Void)
import Numeric.Natural (Natural)

{- | @Injection@ describes a lossless conversion that includes one type in another.

The sole method of this class,

> inject :: from -> into

takes a value @input :: from@ and returns a value @output :: into@ which preserves all the information contained in the input.
Specifically, each @input@ is mapped to a /unique/ @output@.
In mathematical terminology, @inject@ is /injective/:

> inject a ≡ inject b → a ≡ b

The name of the class is derived from the mathematical term.

@Injection@ models the "is-a" relationship used in languages with subtypes (such as in object-oriented programming),
but an explicit cast with @inject@ is required in Haskell.

Although it is often possible to infer the type parameters of this class,
it is advisable to specify one or both of the parameters to @inject@
using a type signature or the @TypeApplications@ language extension.
Specifying the type parameters will give clearer error messages from the type checker in any case.

-}
class Injection from into where
    inject :: from -> into

{- | @Retraction@ undoes an 'Injection'.

Because 'Injection' is a lossless conversion, we can define a @Retraction@ which undoes it.
The method

> retract :: into -> Maybe from

is the (left) inverse of 'inject':

> retract (inject x) = Just x

'retract' is partial (returns 'Maybe') because the type @into@ may be larger than the type @from@;
that is, there may be values in @into@ which are not 'inject'-ed from @from@,
and in that case @retract@ may return 'Nothing'.

Although it is often possible to infer the type parameters of this class,
it is advisable to specify one or both of the parameters to @retract@
using a type signature or the @TypeApplications@ language extension.
Specifying the type parameters will give clearer error messages from the type checker in any case.

-}
class Injection from into => Retraction from into where
    retract :: into -> Maybe from

instance Injection a a where
    inject = id
    {-# INLINE inject #-}

instance Retraction a a where
    retract = Just
    {-# INLINE retract #-}

instance Typeable a => Injection a Dynamic where
    inject = toDyn
    {-# INLINE inject #-}

instance Typeable a => Retraction a Dynamic where
    retract = fromDynamic
    {-# INLINE retract #-}

instance Injection a b => Injection a (Maybe b) where
    inject = Just . inject
    {-# INLINE inject #-}

instance Retraction a b => Retraction a (Maybe b) where
    retract = \x -> x >>= retract @a @b
    {-# INLINE retract #-}

instance Injection a b => Injection (Maybe a) [b] where
    inject = maybeToList . fmap (inject @a @b)
    {-# INLINE inject #-}

instance Retraction a b => Retraction (Maybe a) [b] where
    retract [] = Just Nothing
    retract [b] = Just <$> retract @a @b b
    retract _ = Nothing

instance Injection Natural Integer where
    inject = toInteger
    {-# INLINE inject #-}

instance Retraction Natural Integer where
    retract x
        | x < 0 = Nothing
        | otherwise = Just (fromInteger x)
    {-# INLINE retract #-}

instance Injection Void any where
    inject = \case {}
    {-# INLINE inject #-}

-- | 'Text.unpack' is the canonical injection @'Text' -> 'String'@.
-- There is no injection in the other direction because 'String' can represent
-- invalid surrogate code points, but 'Text' cannot; for details, see
-- "Data.Text".
instance Injection Text String where
    inject = Text.unpack
    {-# INLINE inject #-}

-- | 'Text.Lazy.unpack' is the canonical injection @'Lazy.Text' -> 'String'@.
-- There is no injection in the other direction because 'String' can represent
-- invalid surrogate code points, but 'Lazy.Text' cannot; for details, see
-- "Data.Text.Lazy".
instance Injection Lazy.Text String where
    inject = Text.Lazy.unpack
    {-# INLINE inject #-}

instance Injection Text Lazy.Text where
    inject = Text.Lazy.fromStrict
    {-# INLINE inject #-}

instance Injection Lazy.Text Text where
    inject = Text.Lazy.toStrict
    {-# INLINE inject #-}

instance HasResolution a => Injection Integer (Fixed a) where
    inject = fromInteger
    {-# INLINE inject #-}

instance HasResolution a => Retraction Integer (Fixed a) where
    retract x = retract @Integer (toRational x)
    {-# INLINE retract #-}

instance Injection a (Const a b) where
    inject = Const
    {-# INLINE inject #-}

instance Injection (Const a b) a where
    inject = getConst
    {-# INLINE inject #-}

instance Injection Integer (Ratio Integer) where
    inject = fromInteger
    {-# INLINE inject #-}

instance Retraction Integer (Ratio Integer) where
    retract x
        | Ratio.denominator x == 1 = Just (Ratio.numerator x)
        | otherwise = Nothing
    {-# INLINE retract #-}

instance Num a => Injection a (Complex a) where
    inject = (:+ 0)
    {-# INLINE inject #-}

instance (Eq a, Num a) => Retraction a (Complex a) where
    retract (x :+ y)
      | y == 0 = Just x
      | otherwise = Nothing
    {-# INLINE retract #-}

instance Injection a (Identity a) where
    inject = Identity
    {-# INLINE inject #-}

instance Injection (Identity a) a where
    inject = runIdentity
    {-# INLINE inject #-}

instance Injection (NonEmpty a) [a] where
    inject (x :| xs) = x : xs
    {-# INLINE inject #-}

instance Retraction (NonEmpty a) [a] where
    retract (x : xs) = Just (x :| xs)
    retract [] = Nothing
    {-# INLINE retract #-}

instance Injection a (Down a) where
    inject = pure
    {-# INLINE inject #-}

instance Injection (Down a) a where
    inject = \(Down a) -> a
    {-# INLINE inject #-}

instance Injection a (Product a) where
    inject = pure
    {-# INLINE inject #-}

instance Injection (Product a) a where
    inject = getProduct
    {-# INLINE inject #-}

instance Injection a (Sum a) where
    inject = pure
    {-# INLINE inject #-}

instance Injection (Sum a) a where
    inject = getSum
    {-# INLINE inject #-}

instance Injection a (Dual a) where
    inject = pure
    {-# INLINE inject #-}

instance Injection (Dual a) a where
    inject = getDual
    {-# INLINE inject #-}

instance Injection a (Monoid.Last a) where
    inject = pure
    {-# INLINE inject #-}

instance Retraction a (Monoid.Last a) where
    retract = Monoid.getLast
    {-# INLINE retract #-}

instance Injection a (Monoid.First a) where
    inject = pure
    {-# INLINE inject #-}

instance Retraction a (Monoid.First a) where
    retract = Monoid.getFirst
    {-# INLINE retract #-}

instance Injection a (Semigroup.First a) where
    inject = pure
    {-# INLINE inject #-}

instance Injection (Semigroup.First a) a where
    inject = Semigroup.getFirst
    {-# INLINE inject #-}

instance Injection a (Semigroup.Last a) where
    inject = pure
    {-# INLINE inject #-}

instance Injection (Semigroup.Last a) a where
    inject = Semigroup.getLast
    {-# INLINE inject #-}

instance Injection a (Max a) where
    inject = pure
    {-# INLINE inject #-}

instance Injection (Max a) a where
    inject = getMax
    {-# INLINE inject #-}

instance Injection a (Min a) where
    inject = pure
    {-# INLINE inject #-}

instance Injection (Min a) a where
    inject = getMin
    {-# INLINE inject #-}

instance Injection a (r -> a) where
    inject = const
    {-# INLINE inject #-}

instance Injection Bool Any where
    inject = Any
    {-# INLINE inject #-}

instance Injection Any Bool where
    inject = getAny
    {-# INLINE inject #-}

instance Injection Bool All where
    inject = All
    {-# INLINE inject #-}

instance Injection All Bool where
    inject = getAll
    {-# INLINE inject #-}

max_stars_repo_name: ttuegel/injection

max_issues_repo_name: ttuegel/injection

max_forks_repo_name: ttuegel/injection


text:
module Main where

import MnistLoader
import Network

import Numeric.LinearAlgebra as L

main :: IO ()
main = testMnist

testMnist = do
  net <- initNet [784, 30, 10]
  exemplars <- loadDataWrapper
  test <- loadTestWrapper

  -- let eta = 0.050
  -- let batchSize = 20
  -- let epochs = 30

  let eta = 0.5
  let batchSize = 16
  let epochs = 30

  putStrLn "Running Stochastic Gradient Descent "
  putStrLn $ "eta: " ++ (show eta)
  putStrLn $ "Batch Size: " ++ (show batchSize)
  putStrLn $ "Epochs: " ++ (show epochs)

  trainedNetwork <- sgd net exemplars test epochs 0 batchSize eta
  return ()

test = do
  net <- initNet [3, 5]

  putStrLn "The Net"
  putStrLn $ show net

  putStrLn "Feed Forward"
  putStrLn $ show $ feedForward net (vect 3 [0, 1, 0])

  putStrLn "A ForwardPass"
  putStrLn $ show $ forwardPass net (vect 3 [0, 1, 0])

  putStrLn "Entire backprop"
  let eb = backprop net (vect 3 [0, 1, 0]) (vect 5 [1, 0, 1, 0, 1])
  putStrLn $ show $ length eb
  putStrLn $ show $ eb

  putStrLn "A Minibatch"
  let x = vect 3 [0,1,0]
  let y = vect 5 [1, 0, 1, 0, 1]
  let batch = [(x, y), (x, y), (x, y)]
  let newn = updateMiniBatch net batch 3.0
  putStrLn $ show newn
  return newn

vect n = (n><1)

max_stars_repo_name: madsbuch/haskell-nn

max_issues_repo_name: madsbuch/haskell-nn

max_forks_repo_name: madsbuch/haskell-nn


text:
-- |
-- Module      : Tiler.hs
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created July 9 2015

-- TODO | -
--        -

-- SPEC | -
--        -

module Tiler (renderPathWithJoints) where

---------------------------------------------------------------------------------------------------
-- We'll need these
---------------------------------------------------------------------------------------------------
import Control.Monad (mapM_, forM_)
import Data.Complex

import qualified Graphics.Rendering.Cairo as Cairo

import qualified Southpaw.Picasso.Palette as Palette
import Southpaw.Utilities.Utilities (pairwise)

---------------------------------------------------------------------------------------------------
-- Data
---------------------------------------------------------------------------------------------------
π = pi
τ = 2*π

---------------------------------------------------------------------------------------------------
-- Types
---------------------------------------------------------------------------------------------------
type Path = [Complex Double]

---------------------------------------------------------------------------------------------------
-- Functions
---------------------------------------------------------------------------------------------------
renderPathWithJoints :: Palette.Colour -> Palette.Colour -> Double -> Double -> Path -> Cairo.Render ()
renderPathWithJoints line joint radius thickness path = do

	-- Render the lines
	Cairo.setLineWidth thickness
	forM_ (pairwise $ path ++ [head path]) $ \ (fr, to) -> do
		Palette.choose line
		Cairo.moveTo (realPart fr) (imagPart fr)
		Cairo.lineTo (realPart to) (imagPart to)
		Cairo.stroke

	-- Render the 'joints'
	forM_ path $ \ (x:+y) -> do
		Palette.choose joint
		Cairo.arc x y radius 0 τ
		Cairo.fill

max_stars_repo_name: SwiftsNamesake/Leopardy

max_issues_repo_name: SwiftsNamesake/Leopardy

max_forks_repo_name: SwiftsNamesake/Leopardy


text:
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE TypeFamilies          #-}
-- |
-- Module     : Data.Matrix.Generic
-- Copyright  : Copyright (c) 2012 Aleksey Khudyakov <alexey.skladnoy@gmail.com>
-- License    : BSD3
-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>
-- Stability  : experimental
--
-- Interface for generic immutable matrices. For matrix transposition
-- and conjugate transposition newtype wrappers are used.
module Data.Matrix.Generic (
    -- * Type class
    IsMatrix(..)
  , Mutable
    -- * Accessors
  , rows
  , cols
  , shape
  , (@!)
  , unsafeIndex
    -- * Converions to/from mutable
  , freeze
  , thaw
  , unsafeFreeze
  , unsafeThaw
    -- * Conversion to string
  , showMatrixWith
    -- * Newtype wrappers
  , Transposed(..)
  , Conjugated(..)
  ) where

import Control.Monad             (liftM)
import Control.Monad.Primitive
import Data.List                 (intercalate)
import Data.Complex              (Complex,conjugate)
import qualified Data.Matrix.Generic.Mutable as M
import           Data.Vector.Generic           (Mutable)

----------------------------------------------------------------
-- Type class
----------------------------------------------------------------

-- | Basic API for immutable matrices. Since there's many way to lay
--   matrix in memory there isn't many operation which work for all of
--   them.
--
--   Methods of this type class shouldn't be used directly.
class M.IsMMatrix (Mutable mat) a => IsMatrix mat a where
  -- | Number of rows
  basicRows :: mat a -> Int
  -- | Number of columns
  basicCols :: mat a -> Int
  -- | Read element from matrix
  basicUnsafeIndex :: mat a -> (Int,Int) -> a
  -- | Convert immutable matrix to mutable. Immutable matrix may not
  --   be used after operation.
  basicUnsafeThaw   :: PrimMonad m => mat a -> m (Mutable mat (PrimState m) a)
  -- | Convert mutable matrix to immutable. Mutable matrix may not be
  --   modified after operation.
  basicUnsafeFreeze :: PrimMonad m => Mutable mat (PrimState m) a -> m (mat a)

----------------------------------------------------------------
-- Accessors
----------------------------------------------------------------

-- | Number of rows
rows :: IsMatrix mat a => mat a -> Int
{-# INLINE rows #-}
rows = basicRows

-- | Number of columns
cols :: IsMatrix mat a => mat a -> Int
{-# INLINE cols #-}
cols = basicCols

-- | Shape of the matrix.
shape :: IsMatrix mat a => mat a -> (Int,Int)
{-# INLINE shape #-}
shape m = (rows m, cols m)

-- | Indexing operator without range checking.
unsafeIndex :: IsMatrix mat a
            => mat a            -- ^ Matrix
            -> (Int,Int)        -- ^ (row,column)
            -> a
{-# INLINE unsafeIndex #-}
unsafeIndex = basicUnsafeIndex

-- | Indexing operator with range checking
(@!) :: IsMatrix mat a
     => mat a                   -- ^ Matrix
     -> (Int,Int)               -- ^ (row,column)
     -> a
{-# INLINE (@!) #-}
m @! a@(i,j)
  | i < 0 || i >= rows m = error "Data.Matrix.Generic.@!: row index out of range"
  | j > 0 || j >= cols m = error "Data.Matrix.Generic.@!: column index out of range"
  | otherwise            = unsafeIndex m a

-- | Convert mutable matrix to immutable.
freeze :: (PrimMonad m, IsMatrix mat a) => Mutable mat (PrimState m) a -> m (mat a)
{-# INLINE freeze #-}
freeze m = do
  unsafeFreeze =<< M.clone m

-- | Convert immutable matrix to mutable.
thaw :: (PrimMonad m, IsMatrix mat a) => mat a -> m (Mutable mat (PrimState m) a)
{-# INLINE thaw #-}
thaw m = do
  M.clone =<< unsafeThaw m

-- | Convert mutable matrix to immutable. Mutable matrix may not be
--   modified after operation.
unsafeFreeze :: (PrimMonad m, IsMatrix mat a) => Mutable mat (PrimState m) a -> m (mat a)
{-# INLINE unsafeFreeze #-}
unsafeFreeze = basicUnsafeFreeze

-- | Convert immutable matrix to mutable. Immutable matrix may not
--   be used after operation.
unsafeThaw :: (PrimMonad m, IsMatrix mat a) => mat a -> m (Mutable mat (PrimState m) a)
{-# INLINE unsafeThaw #-}
unsafeThaw = basicUnsafeThaw

-- | Generic function for printing matrix. Mostly useful for debugging
--   purposes.
showMatrixWith :: IsMatrix m a => (a -> String) -> m a -> String
showMatrixWith f m
  = unlines
  $ (show (rows m) ++ " >< " ++ show (cols m))
  : [ intercalate "\t" [f $ unsafeIndex m (i,j) | j <- [0 .. cols m - 1]]
    | i <- [0 .. rows m - 1]
    ]

----------------------------------------------------------------
-- Newtype wrappers
----------------------------------------------------------------

-- | Transposed matrix or vector. Being newtype this wrapper type is
--   used to select different instances for multiplication.
newtype Transposed mat a = Transposed { unTranspose :: mat a }

type instance Mutable (Transposed mat) = M.TransposedM (Mutable mat)

instance IsMatrix mat a => IsMatrix (Transposed mat) a where
  basicRows (Transposed m) = cols m
  {-# INLINE basicRows #-}
  basicCols (Transposed m) = rows m
  {-# INLINE basicCols #-}
  basicUnsafeIndex (Transposed m) (i,j) = unsafeIndex m (j,i)
  {-# INLINE basicUnsafeIndex #-}
  basicUnsafeThaw (Transposed m) = M.TransposedM `liftM` unsafeThaw m
  {-# INLINE basicUnsafeThaw #-}
  basicUnsafeFreeze (M.TransposedM m) = Transposed `liftM` unsafeFreeze m
  {-# INLINE basicUnsafeFreeze #-}

-- | Conjugate-transposed matrix or vector. Being newtype this wrapper type is
--   used to select different instances for multiplication.
newtype Conjugated mat a = Conjugated { unConjugate :: mat a }

type instance Mutable (Conjugated mat) = M.ConjugatedM (Mutable mat)

instance (IsMatrix mat (Complex a), RealFloat a) => IsMatrix (Conjugated mat) (Complex a) where
  basicRows (Conjugated m) = cols m
  {-# INLINE basicRows #-}
  basicCols (Conjugated m) = rows m
  {-# INLINE basicCols #-}
  basicUnsafeIndex (Conjugated m) (i,j)
    | j >= i    = unsafeIndex m (j,i)
    | otherwise = conjugate $! unsafeIndex m (j,i)
  {-# INLINE basicUnsafeIndex #-}
  basicUnsafeThaw (Conjugated m) = M.ConjugatedM `liftM` unsafeThaw m
  {-# INLINE basicUnsafeThaw #-}
  basicUnsafeFreeze (M.ConjugatedM m) = Conjugated `liftM` unsafeFreeze m
  {-# INLINE basicUnsafeFreeze #-}

max_stars_repo_name: Shimuuar/blas-lapack

max_issues_repo_name: Shimuuar/blas-lapack

max_forks_repo_name: Shimuuar/blas-lapack


text:
-- | Now we plot!

module QuickHullDraw where

import           Control.Arrow
import           Data.List
import qualified Data.Vector             as V
import qualified Data.Vector.Unboxed     as UV
import           Graphics.Rendering.Plot
import           Numeric.LinearAlgebra   (Vector (..))
import           QuickHull
import           System.Environment
import           System.Random

main :: IO ()
main = do
  [d] <- getArgs
  let vectors =
        V.fromList .
        map UV.fromList .
        chop 2 . take ((* 2) $ read d :: Int) . randomRs (0 :: Double, 100) $
          mkStdGen 0
      (rxL, rxH) = (0, 100)
      (ryL, ryH) = (0, 100)
      (vx, vy) = convertPoints vectors
      hull' = qhull2D vectors
      (hx, hy) =
        convertPoints . V.fromList $ hull'
  writeFigure PNG ("QHullTest" ++ d ++ "Points" ++ ".png") (600, 600) (testFig vx vy hx hy (rxL, rxH) (ryL, ryH))

convertPoints :: V.Vector (UV.Vector Double) -> (Vector Double, Vector Double)
convertPoints =
  (V.convert *** V.convert) .
  V.unzip . V.map (UV.unsafeHead &&& UV.unsafeLast)

testFig :: Series
        -> Series
        -> Series
        -> Series
        -> (Double, Double)
        -> (Double, Double)
        -> Figure ()
testFig xs ys hx hy (rxL, rxH) (ryL, ryH) = do
  withTextDefaults $ setFontFamily "OpenSymbol"
  withTitle $ setText "2D Convex Hull with QuickHull!"
  setPlots 1 1
  withPlot (1, 1) $
    do
      setDataset [(xs, point ys (Asterisk, red)), (hx, point hy (Box, purple))]
      addAxis XAxis (Side Lower) $ withAxisLabel $ setText "x-axis"
      addAxis YAxis (Side Lower) $ withAxisLabel $ setText "y-axis"
      addAxis XAxis (Value 0) $ return ()
      setRange XAxis Lower Linear (rxL - 5) (rxH + 5)
      setRange YAxis Lower Linear (ryL - 5) (ryH + 5)

chop :: Int -> [a] -> [[a]]
chop n =
  unfoldr
    (\v -> if null v
             then Nothing
             else Just $ splitAt n v)

max_stars_repo_name: emmanueldenloye/QuickHull

max_issues_repo_name: emmanueldenloye/QuickHull

max_forks_repo_name: emmanueldenloye/QuickHull


text:
----------------------
-- 2018.1
-- sule
-- mnist neutral network for AI course homework
----------------------
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}

module Main where

import Data.List(foldl')
import qualified Data.Array as A
import Control.Monad
import qualified Data.ByteString.Char8 as B
import qualified Data.ByteString.Lazy.Char8 as L
import Data.Serialize
import qualified Data.Vector.Storable as V
import Grenade
import Numeric.LinearAlgebra (maxIndex)
import qualified Numeric.LinearAlgebra.Static as S
import System.Random
import System.Directory
import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Diagrams(toFile)
import qualified Data.Array.IO as MA

-- |type of input and output of the network
type IOShape = (S ('D2 28 28), S ('D1 10))

-- |type of the neutral network
type HandRC
  = Network
  '[ Convolution 1 12 5 5 1 1, Pooling 2 2 2 2, Relu
     , Convolution 12 24 5 5 1 1, Pooling 2 2 2 2
     , Reshape, Relu
     , FullyConnected 384 80, Logit
     , FullyConnected 80 10, Logit]
  '[ 'D2 28 28, 'D3 24 24 12, 'D3 12 12 12, 'D3 12 12 12
     , 'D3 8 8 24, 'D3 4 4 24, 'D1 384, 'D1 384
     , 'D1 80, 'D1 80, 'D1 10, 'D1 10]

-- |parse the input and label from a csv file
parseData :: FilePath -> IO [IOShape]
parseData path = do
  strs <- tail . L.lines <$!> L.readFile path
  return (map proc strs)
  where
    proc l =
      ( fromJust $! fromStorable $! V.fromList $! tail lbs
      , fromJust $! fromStorable $! V.fromList $! conv $! head lbs)
      where
        lbs = map (read . L.unpack) $! L.split ',' l
    fromJust (Just a) = a
    conv x' =
      let x = floor x'
      in replicate x 0 ++ [1] ++ replicate (9 - x) 0

-- |train the network with a pair of input and label
trainOne :: LearningParameters -> HandRC -> IOShape -> HandRC
trainOne rate !net p@(d, _) = uncurry (train rate net) p

-- |run the network, get the predict
run :: HandRC -> S ('D2 28 28) -> S ('D1 10)
run = runNet

-- |helper function for random sort
randomSortT :: [Int] -> IO [Int]
randomSortT [] = return []
randomSortT !ls = do
  a <- randomRIO (0, length ls - 1)
  (\ls' -> (ls !! a) : ls') <$!> (randomSortT $! (take a ls ++ drop (a + 1) ls))

-- |helper function for random sort
concatBy::[Int]->[[a]]->[a]
concatBy [] _ = []
concatBy (x:xs) as = (as!!x)++concatBy xs as

-- |implement of random sort
randomSort :: [Int]->Int -> IO [Int]
randomSort ls step=
  if  lg<= step
    then randomSortT ls
    else do
      let ts = lg `div` step
      xr<-randomSortT [0..ts-1]
      res''<-mapM (randomSortT . get) [1 .. div (length ls) step]
      let res'=concatBy xr res''
      if lg `mod` step==0
        then return res'
        else do
          t<-randomRIO (0,ts)
          res<-randomSort (drop (ts*step) ls) step
          return $! take t res'++ res ++drop t res'
  where
    get x = drop (step * (x - 1)) $! take (step * x) ls
    lg =length ls

-- |train one generation
trainOnce ::
     Int->LearningParameters -> Int -> [IOShape] -> [Int] -> HandRC -> IO HandRC
trainOnce gx lp time dts ls !net = do
  let lgT = length ls
  print $! lgT
  foldM
    (\n i -> do
        let ir = ls !! i
        print $! "training the " ++ show i ++ " th data " ++ show ir
        when (i `mod` gx==0) $! do
          let trained = take i dts
          let testD = drop lgT dts
          print $! "testing " ++ show i ++ " records trained before"
          let trainCR = validateList n trained
          appendFile "trained" $! Prelude.show (i, correctRate trainCR)++"\n"
          print $! "testing " ++ show (length dts-lgT) ++ " records for validate"
          let trainTest = validateList n testD
          appendFile "test" $! Prelude.show (i, correctRate trainTest)++"\n"
          plotFromFile ("trained", "test")
        return $!
         trainOne (lp {learningRate = learningRate lp * 0.7 ^ time}) n (dts !! ir))
    net
    [0 .. lgT - 1]

-- |train the nerwork with the given init rate, generations, time to calculate accuracy, csv path, and proportion to train
-- example: trainT 0.007 15 1000 "train.csv" 0.8
-- means train the "train.csv" with a init rate of 0.007, keep training for 15 generations, per 1000 records calculate accuracy
-- if claim a pre-trained network in the directory, it will be loaded
trainT :: FilePath->Double->Int ->Int-> FilePath -> Double -> IO ()
trainT preNet initRate time gx path trainS = do
  dts <- parseData path
  let lg = length dts - 1
  let lgT = 1 + floor (trainS * fromIntegral lg)
  bl<-doesFileExist preNet
  net <- if bl
          then print"loaded">>load preNet
          else randomNetwork
  let lp = LearningParameters initRate 0.9 0.0005
  net' <-
    foldM
      (\nt t -> do
         lsT <- randomSort [0 .. lgT-1] 5
         lsV <- randomSort [lgT .. lg] 5
         trainOnce gx lp t dts lsT nt)
      net
      [1 .. time]
  save net' "net"

calcConfuseMatrix::IO ()
calcConfuseMatrix = do
  net<-load "net"
  dts<-parseData "train.csv"
  ar<-MA.freeze=<<confuseMatrix net dts::IO (A.Array (Int,Int) Int)
  writeFile "confuseMatrix" $ renderA ar

renderA::A.Array (Int,Int) Int->String
renderA = foldl' (\str ((_,a),e)->if a==0 then str++"\n"++show e else str++" "++show e) [] . A.assocs

-- |calculate the confuseMatrix using the given network and records
confuseMatrix::HandRC->[IOShape]->IO (MA.IOUArray (Int,Int) Int)
confuseMatrix net ls = do
  ar<-MA.newArray ((0,0),(9,9)) 0
  mapM_ (\((dt,lb),i)->do
    print i
    let i=(getIdx lb, getIdx $run net dt)
    x<-MA.readArray ar i
    MA.writeArray ar i $x+1) $zip ls [1..]
  return ar

-- |k-cross validation
kCross::Int->Int->Double->FilePath->IO ()
kCross gen k initRate path = do
  lds<-parseData path
  let ws = length lds `div` k
  mapM_ (cross lds ws) [0..k-1]
  where cross lds ws t= do
        let w=t*ws
        let trainD = take w lds++drop (w+ws) lds
        let valD = drop w $! take (w+ws) lds
        net<-trainList initRate t gen trainD
        let tr=validateList net trainD
        let vr=validateList net valD
        save net ("k-cross-net-"++show t)
        putStrLn $!"Train "++show (correctRate tr)++" Validation "++show (correctRate vr)
        writeFile ("k-cross-accuracy-"++show t) $!"Train "++show (correctRate tr)++" Validation "++show (correctRate vr)

-- |helper function for k-cross
justTrain :: LearningParameters -> Int -> [IOShape] -> HandRC -> IO HandRC
justTrain lp time dts !net = do
  let lgT = length dts
  foldM
    (\n i -> do
        print $! "training the " ++ show i ++ " th data"
        return $!
         trainOne (lp {learningRate = learningRate lp * 0.9 ^ time}) n (dts !! i))
    net
    [0 .. lgT - 1]

-- |train a list of records, helper function for k-cross
trainList :: Double->Int->Int->[IOShape] -> IO HandRC
trainList initRate r gen dts = do
  let lg = length dts
  print lg
  net <- randomNetwork :: IO HandRC
  let lp = LearningParameters initRate 0.9 0.0005
  foldM
      (\nt t -> do
         print $! "generation "++show t
         justTrain lp t dts nt)
      net
      [1 .. gen]

-- |get the correctRate from a list of bool
correctRate :: [Bool] -> Double
correctRate [] = 0
correctRate !bs =
  fromIntegral (length (filter id bs)) / fromIntegral (length bs)

-- |validate a list of records
validateList :: HandRC -> [IOShape] -> [Bool]
validateList _ [] = []
validateList !net !ls = map (validateOne net) ls

-- |validate one records
validateOne :: HandRC -> IOShape -> Bool
validateOne !net (!dt, !lb) = judge (run net dt) lb

validate :: IOShape -> IOShape -> Bool
validate (_, !lb') (_, !lb) = judge lb' lb

judge :: S ('D1 10) -> S ('D1 10) -> Bool
judge !x !y = getIdx x == getIdx y

getIdx :: S ('D1 10) -> Int
getIdx (S1D !a) = maxIndex (S.extract a)

-- |save a network to file
save :: HandRC -> FilePath -> IO ()
save net path = B.writeFile path $! runPut (put net)

load :: FilePath -> IO HandRC
load path = do
  modelData <- B.readFile path
  either fail return $! runGet (get :: Get HandRC) modelData

plotFromFile::(FilePath,FilePath)->IO ()
plotFromFile (t,v)=do
  ts<-map read.lines<$>readFile t::IO [(Int,Double)]
  vs<-map read.lines<$>readFile t::IO [(Int,Double)]
  plotAccuracy (trans ts) "train accuracy" "train-accuracy.svg"
  plotAccuracy (trans vs) "validate accuracy" "validate-accuracy.svg"
  where trans xs = filter (\(a,b)->b/=0) $map (\((a,b),z)->(a+z,b)) $ zip xs $concat [replicate 34 x| x<-[0,34000..]]

plotAccuracy::[(Int,Double)]->String->FilePath->IO ()
plotAccuracy xs am path = toFile def path $ do
    layout_title .= "Accuracy"
    setColors [opaque green]
    plot (line am [xs])

runT::IO ()
runT = do
  putStrLn "Enter the pre-trained neutral network file(if not exist then init randomly)"
  net<-getLine
  putStrLn "Enter the init learning rate"
  initRate<-read<$>getLine
  putStrLn "Enter the generations to train(how many repetition)"
  time<-read<$>getLine
  putStrLn "Enter the duration of validation"
  gx<-read<$>getLine
  putStrLn "Enter the percentage for training(1-p for validation)"
  trainS<-read<$>getLine
  putStrLn "Enter the path of training data(train.csv)"
  path<-getLine
  trainT net initRate time gx path trainS

runK::IO ()
runK = do
  putStrLn "Enter the k for k-cross validation"
  k<-read<$>getLine
  putStrLn "Enter the init learning rate"
  initRate<-read<$>getLine
  putStrLn "Enter the generations to train(how many repetition)"
  time<-read<$>getLine
  putStrLn "Enter the path of training data(train.csv)"
  path<-getLine
  kCross time k initRate path

main = do
  putStrLn "choose one action"
  putStrLn "0 -- run randomly training and validation"
  putStrLn "1 -- run k cross validation"
  x<-read<$>getLine::IO Int
  if x==0
    then runT
    else if x==1
      then runK
      else putStrLn "wrong option">>main


text:
module Model where

import Control.Monad.State.Strict
--import Control.Monad.State
import System.Random (StdGen, mkStdGen, randomRs)
import Numeric.LinearAlgebra (Matrix(..), (><), toLists)

type Ising = State IsingState
type Spin = Double
type Model = Matrix Spin

data IsingState = IsingState
     { dim :: Int
     , j :: Double
     , h :: Double
     , step :: Int
     , nAccept :: Int
     , propFreq :: Int
     , model :: Model
     , rng :: StdGen
     }

instance Show IsingState where
    show state =
        "IsingState { d = " ++ show (dim state)
            ++ ", j = " ++ show (j state)
            ++ ", step = " ++ show (step state)
            ++ ", accept = " ++ show (nAccept state)
            ++ ", rng = " ++ show (rng state)
            ++ " }"
            ++ showModel (model state)

showModel :: Model -> String
showModel m =
    let rows = toLists m
        showRow r = "|" ++ foldl (\acc s -> acc ++ showSpin s) "" r ++ "|"
        showSpin s = if s == -1 then " " else "*"
        stringRows = map showRow rows
    in foldl (\acc r -> acc ++ "\n" ++ r) "" stringRows

newModel :: Int -> Int -> Double -> [Int] -> IsingState
newModel seed n j spins =
    let rng = mkStdGen seed
        model = (n><n) $ map fromIntegral spins
    in IsingState
        { dim = n
        , j = j
        , h = 0.0
        , step = 0
        , nAccept = 0
        , propFreq = 1000
        , model = model
        , rng = rng
        }

randModel seed n j =
    let r = mkStdGen seed
        spins = take (n*n) $ randomRs (0,1) r :: [Int]
    in newModel seed n j [if s == 0 then -1 else s | s <- spins]

upSpins   n = [ 1 | _ <- [1..n*n]]
downSpins n = [-1 | _ <- [1..n*n]]

setPropertyFreq :: Int -> IsingState -> IsingState
setPropertyFreq f model = model { propFreq = f }

setH :: Double -> IsingState -> IsingState
setH h model = model { h = h }