text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE NamedFieldPuns #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}

module Statistics.Array
  ( median
  , minimum
  , maximum
  , quartiles
  , listDerivative
  , mad
  , bowleySkew
  ) where

import Prelude hiding (minimum,maximum)

import Control.Monad (forM_)
import Data.Int (Int64)
import Data.Primitive (Prim,PrimArray)
import Data.Primitive.Contiguous (index, write)
import Statistics.Array.Types (Quartiles(..),AscList(..))

import qualified Data.Primitive.Contiguous as Arr
import qualified Data.Primitive.Sort as Arr
import qualified Statistics.Array.Types as Asc

-- | @O(1)@ Median element.
median :: Prim a
  => AscList a -> a
{-# inlineable median #-}
{-# specialize median :: AscList Int -> Int #-}
{-# specialize median :: AscList Int64 -> Int64 #-}
{-# specialize median :: AscList Double -> Double #-}
median (AscList arr) = Arr.index arr (Arr.size arr `div` 2)

-- | @O(1)@ Minimum element.
minimum :: Prim a
  => AscList a -> a
{-# inlineable minimum #-}
{-# specialize minimum :: AscList Int -> Int #-}
{-# specialize minimum :: AscList Int64 -> Int64 #-}
{-# specialize minimum :: AscList Double -> Double #-}
minimum (AscList arr) = Arr.index arr 0

-- | @O(1)@ Maximum element.
maximum :: Prim a
  => AscList a -> a
{-# inlineable maximum #-}
{-# specialize maximum :: AscList Int -> Int #-}
{-# specialize maximum :: AscList Int64 -> Int64 #-}
{-# specialize maximum :: AscList Double -> Double #-}
maximum (AscList arr) = Arr.index arr (Arr.size arr - 1)

quartiles :: Prim a
  => AscList a -> Quartiles a
quartiles (AscList arr) = Quartiles
  { q1 = Arr.index arr (len `div` 4)
  , q2 = Arr.index arr (len `div` 2)
  , q3 = Arr.index arr (3 * len `div` 4)
  }
  where
  len = Arr.size arr

-- hmmm, I guess I'm going with Int64 here to match with Chronos
-- times will be stored in number of nanoseconds after unix epoch
-- which allows me times up to about the year 2262

-- so:
--   check skew of time diffs
--   check dispersion of time diffs (median absolute deviation)
--   check skew of data sizes
--   check dispersion of time diffs (median absolute deviation)
--   check size of data (at the mode)

-- | Median of absolute deviations.
-- This is a measure of dispersion that is more robust than standard deviation.
mad :: forall a.
     (Prim a, Ord a, Num a)
  => AscList a -> a
{-# inlineable mad #-}
{-# specialize mad :: AscList Int -> Int #-}
{-# specialize mad :: AscList Int64 -> Int64 #-}
{-# specialize mad :: AscList Double -> Double #-}
mad asc@(AscList arr) =
  let m = median asc
      devs = (\x -> abs (x - m)) `Arr.map` arr :: PrimArray a
      devsAsc = Arr.sort devs
   in median $ Asc.unsafeFromAscendingArray devsAsc

-- | Take the difference between adjacent elements. Resulting array has
-- length one less than the argument array.
listDerivative :: (Prim a, Num a)
  => AscList a -> PrimArray a
{-# inlineable listDerivative #-}
{-# specialize listDerivative :: AscList Int -> PrimArray Int #-}
{-# specialize listDerivative :: AscList Int64 -> PrimArray Int64 #-}
{-# specialize listDerivative :: AscList Double -> PrimArray Double #-}
listDerivative (AscList arr) = Arr.create $ do
  let len' = Arr.size arr - 1
  diffs <- Arr.new len'
  let go !dstIx !prevElement = if dstIx < len'
        then do
          let !currentElement = index arr (dstIx + 1)
          write diffs dstIx (currentElement - prevElement)
          go (dstIx + 1) currentElement
        else pure diffs
  go 0 (index arr 0)

bowleySkew :: (Integral a) => Quartiles a -> Double
{-# inline bowleySkew #-}
bowleySkew Quartiles{q1,q2,q3}
  -- the limit as q1 approaches q3 is of course dependant on the direction of approach in the 3d space {q1,q2,q3}
  -- realistically, I'm going to call this zero, since real data has no approach, and the dirac distribution has zero skew
  | q1 == q3 = 0.0
  -- if `q1 == q2`, then skew simplifies to `(q3 - q1) / (q3 - q1)`, or `1`
  | q2 == q1 = 1.0
  -- similarly, if `q2 == q3`, then skew simplifies to `-1`
  | q2 == q3 = -1.0
  -- I'm not sure if either of those situations are reliable indicators of skew, though
  | otherwise = (realToFrac $ q3 + q1 - 2*q2) / (realToFrac $ q3 - q1)

max_stars_repo_name: Zankoku-Okuno/array-statistics

max_issues_repo_name: Zankoku-Okuno/array-statistics

max_forks_repo_name: Zankoku-Okuno/array-statistics


text:
-- file:Astro/AccretionDisk.hs
-- Numerical tools for accretion disk analysis

module Astro.AccretionDisk
	where

import Control.Monad
import System.IO.Unsafe
import System.Random
import Numeric.LinearAlgebra

import Astro.Misc

seed = 1337

rndDoubles :: Int -> (Double,Double) -> IO [Double]
rndDoubles n (min,max) = replicateM n $ randomRIO (min,max)

rndDouble :: (Double, Double) -> IO Double
rndDouble (min, max) = randomRIO (min,max) :: IO Double

-- Random position and corresponding velocity in a flat xy-plane Keplerian disk
-- with CCW rotation when viewed towards negative z-axis
rndDiskPosVel :: (Double, Double) -> Double -> IO (Vector Double, Vector Double)
rndDiskPosVel (minR,maxR) gpar =
	do
	r <- rndDouble (minR,maxR)
	th <- rndDouble (0, 2*pi)
	let rv = fromList [ r * cos th, r * sin th, 0]
	    vv = fromList [-rv@>1, rv@>0, 0.0]
	    vv' = scale (sqrt gpar/norm rv) vv
	return (rv, vv')

-- Create a flat disk (planar), rotating in a CCW sense when viewed
-- towards negative z-axis
flatDisk :: (Double, Double) -> Double -> Int
		-> [(Vector Double, Vector Double)]
flatDisk (minR, maxR) gpar nparts = unsafePerformIO $
	do
		setStdGen $ mkStdGen seed
		replicateM nparts $ rndDiskPosVel (minR,maxR) gpar

ext: hs

max_stars_repo_name: sageh/AstroHaskell

max_issues_repo_name: sageh/AstroHaskell

max_forks_repo_name: sageh/AstroHaskell


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

-- | This module provides the Jordan-Wigner transformation and
-- symbolic derivation of circuit templates for second quantized
-- interaction terms. It is essentially a fully automated version of
-- the calculations from
--
-- * James D. Whitfield, Jacob Biamonte, and Alán
-- Aspuru-Guzik. \"Simulation of electronic structure Hamiltonians
-- using quantum computers.\"
-- /Molecular Physics/ 109(5):735–750, 2011.
-- See also <http://arxiv.org/abs/1001.3855v3>.

module Algorithms.GSE.JordanWigner where

import Quipper

import QuipperLib.Decompose

import Data.Complex
import qualified Data.Map as Map
import Data.Map (Map)

import Libraries.Auxiliary (sequence_right_)

import Text.Printf

-- ----------------------------------------------------------------------
-- * Overview

-- $ For a given tuple of orbital indices, (/p/,/q/) in case of
-- one-electron interactions, or (/p/,/q/,/r/,/s/) in case of two-electron
-- interactions, we first calculate the Jordan-Wigner transformation
-- of the second quantized hermitian interaction terms
--
-- /a/[sub /p/][sup †]/a/[sub /p/],
--
-- /a/[sub /p/][sup †]/a/[sub /q/] + /a/[sub /q/][sup †]/a/[sub /p/],
--
-- /a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /q/]/a/[sub /p/],
--
-- /a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /r/]/a/[sub /s/] +
-- /a/[sub /s/][sup †]/a/[sub /r/][sup †]/a/[sub /q/]/a/[sub /p/].
--
-- Next, we decompose each operator into a linear combination /H/ =
-- λ[sub 1]/M/[sub 1] + ... + λ[sub /n/]/M/[sub /n/] of mutually
-- commuting hermitian tensors. At this point, each summand /M/[sub /j/]
-- in the linear combination will be a tensor product of the following
-- operators (not necessarily in this order):
--
-- * an even number (possibly zero) of Pauli /X/ operators;
--
-- * an even number (possibly zero) of Pauli /Y/ operators;
--
-- * zero or more Pauli /Z/ operators, and
--
-- * zero or more /D/ operators, where /D/ = σ[sup −]σ[sup +] = (/I/−/Z/)\/2.
--
-- Note that there may be zero terms in the summation; this happens,
-- for example, for
-- /a/[sub /p/][sup †]/a/[sub /p/][sup †]/a/[sub /r/]/a/[sub /s/],
-- because two electrons cannot occupy the same spin orbital due to
-- their fermionic nature. In this case, /H/ = 0 and [exp -/i/θ/H/] = /I/.
--
-- Next, we calculate [exp -/i/θ/H/]. Because the summands /M/[sub /j/]
-- commute, we can exponentiate each summand separately, using the formula
-- [exp -/i/θ/H/] = [exp -/i/θλ[sub 1]/M/[sub 1]]⋯[exp -/i/θλ[sub /n/]/M/[sub /n/]].
--
-- We then generate the circuit for [exp -/i/θλ[sub /j/]/M/[sub /j/]]
-- by applying a sequence of basis changes until the problem is
-- reduced to a controlled rotation. The basis changes are, in this
-- order:
--
-- 1. Change each Pauli /X/ operator in /M/[sub /j/] to a Pauli /Z/
-- operator, and apply a Hadamard basis change to the corresponding
-- qubit. This uses the relation /HXH/ = /Z/.
--
-- 2. Change each Pauli /Y/ operator in /M/[sub /j/] to a Pauli /Z/
-- operator, and apply a [bold Y] basis change to the corresponding
-- qubit. Note: the [bold Y] basis change operator is defined in
-- [Whitfield et al.] as /R/[sub /x/](-π\/2) = (/I/+/iX/)\/√2, or
-- equivalently [bold Y] = /SHS/, and satisfies [bold Y][super †]/Y/[bold Y] =
-- /Z/. It should not be confused with the Pauli /Y/ operator.
--
-- 3. If the operator /M/[sub /j/] contains one or more Pauli /Z/ operators
-- (including those obtained in steps 1 and 2), then do a basis change
-- by a cascade of controlled-not gates to reduce this to a single /Z/
-- operator. This uses the relation /CNot/ (/Z/⊗/Z/) /CNot/ = /I/⊗/Z/.
--
-- After these basis changes, the operator /M/[sub /j/] consists of
-- exactly zero or one Pauli /Z/ operator, together with zero or more
-- /D/ operators.  To see how to translate this into a controlled
-- rotation, note that for any operator /A/, we have
--
-- \[image expDA.png]
--
-- Therefore, each /D/ operator in /M/[sub /j/] turns into a control
-- after exponentiation. The final rotation is then computed by
-- distinguishing two cases:
--
-- * If /M/[sub /j/] contains a Pauli /Z/ operator, then use the
-- relation [exp -/i/θ/Z/] = /R/[sub /z/](2θ). In this case, the circuit
-- for [exp -/i/θ/M/[sub /j/]] is a controlled /R/[sub /z/](2θ) gate in
-- the position of the /Z/ operator, with zero or more controls in the
-- positions of any /D/ operators.
--
-- * If /M/[sub /j/] does not contain a Pauli /Z/ operator, then the
-- operation to be performed is a phase change [exp -/i/θ], controlled
-- by the qubits in the positions of the /D/ operators. Note that
-- there must be at least one /D/ operator in this case. Also note
-- that a controlled [exp -/i/θ] gate is identical to a /T/(θ) gate.

-- ----------------------------------------------------------------------
-- * Correctness of the templates

-- $ As outlined above, the functions in this module generate each
-- circuit from first principles, based on the Jordan-Wigner
-- representation of operators and on algebraic transformations. They
-- do not rely on pre-fabricated circuit templates.
--
-- Based on the automated calculations provided by this module, we
-- have found small typos in the 5 templates provided by [Whitfield
-- et al.] (Table 3, or Table A1 in the arXiv version).
--
-- * The template for the number-excitation operator is missing a
-- control on its rotation gate.
--
-- * In the template for the Coulomb operator, the angles are wrong.
-- Moreover, this program finds a simpler template.
--
-- * In the template for the double excitation operator, the angles
-- are wrong; they should be ±θ\/4 instead of θ.
--
-- The corrected templates generated by our code are as follows:
--
-- * Number operator /h/[sub /pp/] /a/[sub /p/][sup †]/a/[sub /p/].
--
-- > [image b0-template.png]
--
-- * Excitation operator /h/[sub /pq/] /a/[sub /p/][sup †]/a/[sub /q/].
--
-- > [image b1-template.png]
--
-- * Coulomb and exchange operators /h/[sub /pqqp/]
-- /a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /q/]/a/[sub /p/].
--
-- > [image b2-template.png]
--
-- * Number-excitation operator /h/[sub /pqqr/]
-- (/a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /q/]/a/[sub /r/] +
-- /a/[sub /r/][sup †]/a/[sub /q/][sup †]/a/[sub /q/]/a/[sub /p/]).
-- The sign of ±θ depends on the relative ordering of the indices /p,q,r/.
--
-- > [image b3-template.png]
--
-- * Double excitation operator
-- /h/[sub /pqrs/]
-- (/a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /r/]/a/[sub /s/] +
-- /a/[sub /s/][sup †]/a/[sub /r/][sup †]/a/[sub /q/]/a/[sub /p/]).
-- The sign of ±θ\/4 in each of the eight terms depends on the relative ordering of the indices /p,q,r,s/.
--
-- > [image b4-template.png]

-- ----------------------------------------------------------------------
-- * Alternate Coulomb templates

-- $ As noted above, our algorithm found the following template for the
-- Coulomb operator
-- /a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /q/]/a/[sub /p/]:
--
-- > [image b2-template.png]
--
-- This is simpler than the template given in [Whitfield et al.], even
-- after one accounts for the cost of decomposing the additional
-- controlled /T/(θ) gate into elementary gates. However, an
-- equivalent circuit can also be given that is more similar to the
-- one in [Whitfield et al.] (but with corrected rotation angles):
--
-- > [image b2-orthodox.png]
--
-- We call this the \"orthodox\" template, because it is closer to the
-- one specified by Whitfield et al. The program will use the orthodox
-- template if the command line option @--orthodox@ is given, and it
-- will use the simplified template otherwise.

-- ----------------------------------------------------------------------
-- * General-purpose auxiliary functions

-- | Construct a list consisting of /n/ repetitions of some element.
power :: Int -> a -> [a]
power n x = take n $ repeat x

-- | Extract a list of /n/-1 consecutive pairs from an /n/-element list:
--
-- > consecutive_pairs [] = []
-- > consecutive_pairs [1] = []
-- > consecutive_pairs [1,2] = [(1,2)]
-- > consecutive_pairs [1,2,3] = [(1,2),(2,3)]
-- > consecutive_pairs [1,2,3,4] = [(1,2),(2,3),(3,4)]
consecutive_pairs :: [a] -> [(a,a)]
consecutive_pairs [] = []
consecutive_pairs [h] = []
consecutive_pairs (h1:h2:t) = (h1,h2) : consecutive_pairs (h2:t)

-- ----------------------------------------------------------------------
-- * Scalars

-- | The type of complex numbers. Here, we use a floating point
-- representation, although a symbolic representation would also be
-- possible. Since for the purpose of this algorithm, all denominators
-- are powers of 2, the floating point representation is in fact exact.
type Scalar = Complex Double

-- | The complex number /i/.
i :: Scalar
i = 0 :+ 1

-- ----------------------------------------------------------------------
-- * Basic Gates

-- | Apply a /R/[sub z](θ)=[exp -/i/θ/Z/\/2] gate. The parameter θ is a
-- Bloch sphere angle.
--
-- \[image Rz.png]
rotZ_at :: Double -> Qubit -> Circ ()
rotZ_at theta q = named_rotation_at "Rz(%)" theta q

-- | Apply a /G/(θ) gate. This is a global phase change of [exp -/i/θ],
-- so this gate only \"does\" something when it is controlled.
-- Although it is logically a 0-ary gate, we give it a qubit argument
-- to specify where the gate can be drawn in circuit diagrams.
--
-- \[image G.png]
gse_G_at :: Double -> Qubit -> Circ ()
gse_G_at theta q = named_rotation_at "G(%)" theta q

-- | Apply a /T/(θ) gate. This is a /Z/-rotation, but differs
-- from /R/[sub z](-θ) by a global phase.
--
-- \[image T.png]
gse_T_at :: Double -> Qubit -> Circ ()
gse_T_at theta q = named_rotation_at "T(%)" theta q

-- | Apply a [bold Y] basis change gate. This is defined as [bold Y] = /SHS/,
-- or equivalently,
--
-- \[image Y.png]
--
-- This should not be confused with the Pauli /Y/ gate.
gse_Y_at :: Qubit -> Circ ()
gse_Y_at q = named_gate_at "YY" q

-- ----------------------------------------------------------------------
-- * Basic operators

-- | This type provides a symbolic representation of certain
-- operators, generated by the Pauli operators, /P/ = σ[sup +], and
-- /M/ = σ[sup −]. For lack of a better term, we call these the
-- \"basic\" operators. Note that apart from /P/ and /M/, all of these
-- are hermitian.
data Op =
  I -- ^ Identity operator.
  | X  -- ^ Pauli /X/ operator.
  | Y  -- ^ Pauli /Y/ operator.
  | Z  -- ^ Pauli /Z/ operator.
  | P  -- ^ σ[sup +] operator = (0,1;0,0).
  | M  -- ^ σ[sup −] operator = (0,0;1,0).
  | A  -- ^ σ[sup +]σ[sup −] operator = (1,0;0,0).
  | D  -- ^ σ[sup −]σ[sup +] operator = (0,0;0,1).
  deriving (Show, Eq, Ord)

-- | A type to represent scalar multiples. An element of ('Scaled'
-- /a/) is a pair (λ, /x/) of a complex scalar λ and an element /x/ ∈
-- /a/.
data Scaled a = Scaled Scalar a

-- | Multiplication of basic operators. Note that the product of two
-- basic operators is not usually itself a basic operator, but a
-- scalar multiple thereof. This multiplication encodes the algebraic
-- laws of basic operators in symbolic form.

-- Implementation note: the multiplication laws are currently
-- implemented as a long case distinction. Perhaps it could be done
-- more cleverly.
mult :: Op -> Op -> Scaled Op

-- The Pauli group
mult I x = Scaled 1 x
mult x I = Scaled 1 x

mult X X = Scaled 1 I
mult X Y = Scaled i Z
mult X Z = Scaled (-i) Y

mult Y X = Scaled (-i) Z
mult Y Y = Scaled 1 I
mult Y Z = Scaled i X

mult Z X = Scaled i Y
mult Z Y = Scaled (-i) X
mult Z Z = Scaled 1 I

-- P, M, D, and A
mult X P = Scaled 1 D
mult X M = Scaled 1 A
mult X A = Scaled 1 M
mult X D = Scaled 1 P

mult Y P = Scaled i D
mult Y M = Scaled (-i) A
mult Y A = Scaled i M
mult Y D = Scaled (-i) P

mult Z P = Scaled 1 P
mult Z M = Scaled (-1) M
mult Z A = Scaled 1 A
mult Z D = Scaled (-1) D

mult P X = Scaled 1 A
mult M X = Scaled 1 D
mult A X = Scaled 1 P
mult D X = Scaled 1 A

mult P Y = Scaled i A
mult M Y = Scaled (-i) D
mult A Y = Scaled (-i) P
mult D Y = Scaled i M

mult P Z = Scaled (-1) P
mult M Z = Scaled 1 M
mult A Z = Scaled 1 A
mult D Z = Scaled (-1) D

mult P P = Scaled 0 I
mult P M = Scaled 1 A
mult P A = Scaled 0 I
mult P D = Scaled 1 P

mult M P = Scaled 1 D
mult M M = Scaled 0 I
mult M A = Scaled 1 M
mult M D = Scaled 0 I

mult A P = Scaled 1 P
mult A M = Scaled 0 I
mult A A = Scaled 1 A
mult A D = Scaled 0 I

mult D P = Scaled 0 I
mult D M = Scaled 1 M
mult D A = Scaled 0 I
mult D D = Scaled 1 D

-- ----------------------------------------------------------------------
-- * Tensors of basic operators

-- | We use a list of basic operators to represent a tensor
-- product. The convention is that infinitely many identity operators
-- are implicitly appended at the end of the list.
type Tensor = [Op]

-- | Normalize a tensor, by stripping away trailing identities.
normalize_tensor :: Tensor -> Tensor
normalize_tensor [] = []
normalize_tensor (h:t) =
  if h == I && null n
  then []
  else (h:n)
    where n = normalize_tensor t

-- | The identity tensor.
tensor_id :: Tensor
tensor_id = []

-- | Multiply two tensors. This returns a scaled tensor.
mult_tensor :: Tensor -> Tensor -> Scaled Tensor
mult_tensor [] bs = Scaled 1 bs
mult_tensor as [] = Scaled 1 as
mult_tensor (a:as) (b:bs) = Scaled (x*y) (c:cs) where
  Scaled x c = mult a b
  Scaled y cs = mult_tensor as bs

-- | Multiply two scaled tensors.
mult_scaled_tensor :: Scaled Tensor -> Scaled Tensor -> Scaled Tensor
mult_scaled_tensor (Scaled x a) (Scaled y b) = Scaled (x*y*z) c where
  Scaled z c = a `mult_tensor` b

-- ----------------------------------------------------------------------
-- * Linear combinations of tensors

-- | A type to represent complex linear combinations of tensors.
type TensorLC = Map Tensor Scalar

-- | The origin.
lc_zero :: TensorLC
lc_zero = Map.empty

-- | Add a tensor to a linear combination.
lc_insert :: TensorLC -> Scaled Tensor -> TensorLC
lc_insert lc (Scaled lambda t) =
  if newvalue == 0
  then
    Map.delete m lc
  else
    Map.insert m newvalue lc
      where
        m = normalize_tensor t
        newvalue = case Map.lookup m lc of
          Nothing -> lambda
          Just x -> lambda + x

-- | Turn a list of scaled tensors into a 'TensorLC'.
lc_from_list :: [Scaled Tensor] -> TensorLC
lc_from_list = foldl lc_insert lc_zero

-- | Turn a 'TensorLC' into a list of scaled tensors.
lc_to_list :: TensorLC -> [Scaled Tensor]
lc_to_list lc = [Scaled x y | (y,x) <- Map.toList lc]

-- ----------------------------------------------------------------------
-- * Jordan-Wigner representation

-- $ The next two functions provide the Jordan-Wigner representation of
-- (Fock-space) annihilation and creation operators.

-- | Construct the Jordan-Wigner annihilation operator /a/[sub /p/] =
-- /IIIIPZZZZZ.../ for spin-orbital index /p/. The first parameter is
-- /p/, and the second one is /M/ (the number of spin-orbitals).
-- Precondition: 0 ≤ /p/ < /M/.
jw :: Int -> Int -> Scaled Tensor
jw p m = Scaled 1 (power p I ++ [P] ++ power (m-p-1) Z)

-- | Construct the Jordan-Wigner creation operator /a/[sub /p/][sup †]
-- = /IIIIMZZZZ.../ for spin-orbital index /p/.  The first parameter
-- is /p/, and the second one is /M/ (the number of spin-orbitals).
-- Precondition: 0 ≤ /p/ < /M/.
jw_dagger :: Int -> Int -> Scaled Tensor
jw_dagger p m = Scaled 1 (power p I ++ [M] ++ power (m-p-1) Z)

-- ----------------------------------------------------------------------
-- * Second quantized interaction terms

-- ** Simple interaction terms

-- | Construct the one-electron second quantized non-hermitianized
-- interaction term /a/[sub /p/][sup †]/a/[sub /q/].  The parameters
-- are /p,q/.
one_electron_operator_simple :: Int -> Int -> Scaled Tensor
one_electron_operator_simple p q = ap * aq
  where
    ap = jw_dagger p m
    aq = jw q m
    m = maximum [p,q] + 1
    (*) = mult_scaled_tensor

-- | Construct the two-electron second quantized non-hermitianized
-- interaction term
-- /a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /r/]/a/[sub /s/].
-- The parameters are /p,q,r,s/.
two_electron_operator_simple :: Int -> Int -> Int -> Int -> Scaled Tensor
two_electron_operator_simple p q r s = ap * aq * ar * as
  where
    ap = jw_dagger p m
    aq = jw_dagger q m
    ar = jw r m
    as = jw s m
    m = maximum [p,q,r,s] + 1
    (*) = mult_scaled_tensor

-- ** Hermitian interaction terms

-- | Construct
-- /a/[sub /p/][sup †]/a/[sub /q/]
-- if /p/ = /q/, and
-- /a/[sub /p/][sup †]/a/[sub /q/] + /a/[sub /q/][sup †]/a/[sub /p/]
-- otherwise.
one_electron_operator :: Int -> Int -> TensorLC
one_electron_operator p q =
  if p == q
  then lc_from_list [a_pq]
  else lc_from_list [a_pq, a_qp]
    where
      a_pq = one_electron_operator_simple p q
      a_qp = one_electron_operator_simple q p

-- | Construct
-- /a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /r/]/a/[sub /s/]
-- if (/p/,/q/) = (/s/,/r/), and
-- /a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /r/]/a/[sub /s/] +
-- /a/[sub /s/][sup †]/a/[sub /r/][sup †]/a/[sub /q/]/a/[sub /p/]
-- otherwise.
two_electron_operator :: Int -> Int -> Int -> Int -> TensorLC
two_electron_operator p q r s =
  if (p,q) == (s,r)
  then lc_from_list [a_pqrs]
  else lc_from_list [a_pqrs, a_srqp]
    where
    a_pqrs = two_electron_operator_simple p q r s
    a_srqp = two_electron_operator_simple s r q p

-- ----------------------------------------------------------------------
-- * /XYZD/ decomposition

-- | Decompose a basic operator into linear combinations of hermitian basic operators.
-- This uses the relations /P/ = 1\/2 /X/ + /i/\/2 /Y/ and
-- /M/ = 1\/2 /X/ - /i/\/2 /Y/.
decompose_basis :: Op -> [Scaled Op]
decompose_basis P = [Scaled (1/2) X, Scaled (i/2) Y]
decompose_basis M = [Scaled (1/2) X, Scaled (-i/2) Y]
decompose_basis x = [Scaled 1 x]  -- default case: the remaining operators are already hermitian

-- | Decompose a tensor into a linear combination of hermitian
-- tensors. Due to sign alternation, the individual tensors all come
-- out to commute with each other.
decompose_tensor :: Tensor -> TensorLC
decompose_tensor [] = lc_from_list [Scaled 1 tensor_id]
decompose_tensor (h:t) =
  lc_from_list [Scaled (x*y) (g:gs) |
                Scaled x g <- decompose_basis h,
                Scaled y gs <- lc_to_list (decompose_tensor t)]

-- | Decompose a linear combination of tensors into a linear
-- combination of hermitian tensors.
decompose_tensor_lc :: TensorLC -> TensorLC
decompose_tensor_lc lc =
  lc_from_list [ Scaled (x*y) g |
                 Scaled x gs <- lc_to_list lc,
                 Scaled y g <- lc_to_list (decompose_tensor gs)]

-- ----------------------------------------------------------------------
-- * Exponentiation and circuit generation

-- | Given a simple hermitian tensor /H/ and an angle θ, generate a
-- circuit for [exp -/i/θ/H/].  The given list of input qubits is in
-- the same order as the operators in /H/. Precondition: /H/ is made
-- up of zero or more identity operators and one or more of the
-- operators /X/, /Y/, /Z/, and /D/. The last parameter is a list of
-- additional controls.

exponentiate_simple :: Scaled Tensor -> Double -> [Qubit] -> [Qubit] -> Circ ()
exponentiate_simple (Scaled s ms) theta qs ctl = do
  -- First analyze the tensor:
  let
    -- Find all X positions
    xs = [ i | (m, i) <- zip ms [0,1..], m == X ]
    -- Find all Y positions
    ys = [ i | (m, i) <- zip ms [0,1..], m == Y ]
    -- Find all X, Y, Z positions
    zs = [ i | (m, i) <- zip ms [0,1..], m `elem` [X, Y, Z] ]
    -- Find all D positions
    ds = [ i | (m, i) <- zip ms [0,1..], m == D ]
  basischange xs ys zs qs
  rotation alpha ds zs `controlled` ctl
  reverse_generic_imp (basischange xs ys zs) qs
  where
    alpha = theta * realPart s
    basischange :: [Int] -> [Int] -> [Int] -> [Qubit] -> Circ ()
    basischange xs ys zs qs = do
      -- for every X or Y in the operator, apply the appropriate basis
      -- change to change it to Z.
      sequence_ [ hadamard_at (qs !! i) | i <- xs ]
      sequence_ [ gse_Y_at (qs !! i) | i <- ys ]
      -- apply a cascade of c-not operators to all X, Y, or Z in the
      -- operator
      sequence_right_ [ qnot_at (qs !! i0) `controlled` (qs !! i1) | (i0,i1) <- consecutive_pairs zs]
    rotation :: Timestep -> [Int] -> [Int] -> Circ ()
    rotation alpha ds zs = do
      case zs of
        [] -> -- if there are no Z operators, produce a controlled
              -- e^{-iα} gate.
          case ds of
            [] -> error "exponentiate_simple: precondition violated"
            d:ds' -> gse_T_at alpha (qs !! d) `controlled` [qs !! d' | d' <- ds']
        z:zs' -> -- if there are Z operators, produce a controlled
                 -- e^{-iαZ} gate.
          rotZ_at (2*alpha) (qs !! z) `controlled` [qs !! d | d <- ds]

-- | Given a tensor /H/ (already decomposed into commuting simple
-- tensors) and an angle θ, generate a circuit for [exp -/i/θ/H/].
-- The given list of input qubits is in the same order as the
-- operators in /H/.
exponentiate :: TensorLC -> Double -> [Qubit] -> [Qubit] -> Circ ()
exponentiate lc theta qs ctl =
  sequence_ [ exponentiate_simple a theta qs ctl | a <- lc_to_list lc ]

-- ----------------------------------------------------------------------
-- * Generate top-level templates

-- | @'one_electron_circuit' theta p q@: Generate the circuit for the
-- hermitianized one-electron interaction with spin-orbital indices
-- /p/, /q/. More precisely, generate [exp -/i/θ/H/], where
-- /H/ = /a/[sub /p/][sup †]/a/[sub /q/] if /p/ = /q/ and
-- /H/ = /a/[sub /p/][sup †]/a/[sub /q/]
-- + /a/[sub /q/][sup †]/a/[sub /p/] otherwise.
--
-- This function recognizes an important special case: if θ=0.0, don't
-- generate any gates at all. The case θ=0.0 frequently arises because
-- of the conversion from spatial orbitals to spin orbitals.
one_electron_circuit :: Double -> Int -> Int -> [Qubit] -> [Qubit] -> Circ ()
one_electron_circuit 0.0 p q qs ctl = return ()
one_electron_circuit theta p q qs ctl = do
  comment_with_label (printf "ENTER: one_electron_circuit (theta=%.3e, p=%d, q=%d)" theta p q) (qs,ctl) ("qs","ctl")
  exponentiate op theta qs ctl
  comment_with_label "EXIT: one_electron_circuit" (qs,ctl) ("qs","ctl")
  where
    op = decompose_tensor_lc (one_electron_operator p q)

-- | @'two_electron_circuit' theta p q r s@:
-- Generate the circuit for the hermitianized two-electron interaction
-- with spin-orbital indices /p/, /q/, /r/, /s/. More precisely, generate
-- [exp -/i/θ/H/], where
-- /H/ = /a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /r/]/a/[sub /s/] if
-- (/p/,/q/) = (/s/,/r/) and /H/ =
-- /a/[sub /p/][sup †]/a/[sub /q/][sup †]/a/[sub /r/]/a/[sub /s/] +
-- /a/[sub /s/][sup †]/a/[sub /r/][sup †]/a/[sub /q/]/a/[sub /p/]
-- otherwise.
--
-- This function recognizes an important special case: if θ=0.0, don't
-- generate any gates at all. The case θ=0.0 frequently arises because
-- of the conversion from spatial orbitals to spin orbitals.
two_electron_circuit :: Double -> Int -> Int -> Int -> Int -> [Qubit] -> [Qubit] -> Circ ()
two_electron_circuit 0.0 p q r s qs ctl = return ()
two_electron_circuit theta p q r s qs ctl = do
  comment_with_label (printf "ENTER: two_electron_circuit (theta=%.3e, p=%d, q=%d, r=%d, s=%d)" theta p q r s) (qs,ctl) ("qs","ctl")
  exponentiate op theta qs ctl
  comment_with_label "EXIT: two_electron_circuit" (qs,ctl) ("qs","ctl")
  where
    op = decompose_tensor_lc (two_electron_operator p q r s)

-- | Like 'two_electron_circuit', but use the \"orthodox\" circuit
-- template for the Coulomb operator /a/[sub /p/][sup †]/a/[sub
-- /q/][sup †]/a/[sub /q/]/a/[sub /p/]. This generates a circuit using
-- three rotations, similar to [Whitfield et al.], but with corrected
-- angles,
--
-- > [image b2-orthodox.png]
--
-- instead of the simpler circuit that 'two_electron_circuit' would
-- normally generate:
--
-- > [image b2-template.png]

two_electron_circuit_orthodox :: Double -> Int -> Int -> Int -> Int -> [Qubit] -> [Qubit] -> Circ ()
two_electron_circuit_orthodox 0.0 p q r s qs ctl = return ()
two_electron_circuit_orthodox theta p q r s qs ctl | q==r && p==s && p/=q = do
  comment_with_label (printf "ENTER: two_electron_circuit_orthodox(theta=%.3e, p=%d, q=%d, r=%d, s=%d)" theta p q r s) (qs,ctl) ("qs","ctl")
  let pp = qs !! p
      qq = qs !! q
  gse_G_at (theta/4) pp `controlled` ctl
  rotZ_at (-theta/2) pp `controlled` ctl
  rotZ_at (-theta/2) qq `controlled` ctl
  qnot_at qq `controlled` pp  -- control ctl is not needed
  rotZ_at (theta/2) qq `controlled` ctl
  qnot_at qq `controlled` pp  -- control ctl is not needed
  comment_with_label "EXIT: two_electron_circuit_orthodox" (qs,ctl) ("qs","ctl")

-- all other cases aren't Coulomb operators, so fall back to
-- two_electron_circuit.
two_electron_circuit_orthodox theta p q r s qs ctl = do
  two_electron_circuit theta p q r s qs ctl

-- ----------------------------------------------------------------------
-- * Testing

-- $ We provide two functions, accessible via command line options,
-- that allow the user to display individual templates.

-- | Display the circuit for the hermitianized one-electron interaction,
-- with θ=1.
show_one_electron :: Format -> GateBase -> Int -> Int -> IO ()
show_one_electron format gatebase p q =
  print_generic format (decompose_generic gatebase circuit) (replicate (n-m) qubit) where
    circuit qs = one_electron_circuit 1.0 (p-m) (q-m) qs []
    n = maximum [p,q] + 1
    m = minimum [p,q]

-- | Display the circuit for the hermitianized two-electron interaction,
-- with θ=1.
show_two_electron :: Format -> GateBase -> Int -> Int -> Int -> Int -> IO ()
show_two_electron format gatebase p q r s =
  print_generic format (decompose_generic gatebase circuit) (replicate (n-m) qubit) where
    circuit qs = two_electron_circuit 1.0 (p-m) (q-m) (r-m) (s-m) qs []
    n = maximum [p,q,r,s] + 1
    m = minimum [p,q,r,s]

-- | Like 'show_two_electron', but use the \"orthodox\" template for
-- the Coulomb operator.
show_two_electron_orthodox :: Format -> GateBase -> Int -> Int -> Int -> Int -> IO ()
show_two_electron_orthodox format gatebase p q r s =
  print_generic format (decompose_generic gatebase circuit) (replicate (n-m) qubit) where
    circuit qs = two_electron_circuit_orthodox 1.0 (p-m) (q-m) (r-m) (s-m) qs []
    n = maximum [p,q,r,s] + 1
    m = minimum [p,q,r,s]

max_stars_repo_name: fritzo/quipper

max_issues_repo_name: fritzo/quipper

max_forks_repo_name: fritzo/quipper


text:
{-# LANGUAGE ScopedTypeVariables #-}

module Math.Probably.NelderMead where

import Math.Probably.FoldingStats

import Data.Ord
import Data.List
import Data.Maybe

import Debug.Trace

import Foreign.Storable

import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as L

{-mean, m2 :: V.Vector Double
mean = fromList [45.1,10.3]

cov :: Matrix Double
cov = (2><2) $ [5.0, 1,
                1, 1.5]

invCov = inv cov

lndet = log $ det cov

pdf = multiNormalByInv lndet invCov mean

m2 = fromList [55.1,20.3]

main = runRIO $ do
  io $ print $ pdf mean
  ap <- nmAdaMet (defaultAM {nsam =1000000}) pdf (fromList [3.0,0.8])
  {-iniampar <- sample $ initialAdaMet 500 5e-3 (pdf) $ fromList [30.0,8.0]
  io $ print iniampar
  froampar <- runAndDiscard 5000 (show . ampPar) iniampar $ adaMet False (pdf)
  io $ print froampar
  io $ print $ realToFrac (count_accept froampar) / (realToFrac $ count froampar) -}
  let iniSim = genInitial (negate . pdf) 0.1 $ fromList [3.0,0.5]
  io $ mapM_ print iniSim
  let finalSim =  goNm (negate . pdf) 1 iniSim
  io $ print $ finalSim
  io $ print $ hessianFromSimplex (negate . pdf) finalSim -}

type Simplex = [(V.Vector Double, Double)]

{-instance (Storable a, Num a) => Num (V.Vector a) where
   (+) = V.zipWith (+)
   (-) = V.zipWith (-)
   (*) = V.zipWith (*) -}

scale x = V.map (*x)

centroid :: Simplex -> V.Vector Double
centroid points = scale (recip l) $ sum $ map fst points
    where l = fromIntegral $ length points

nmAlpha = 1
nmGamma = 2
nmRho = 0.5
nmSigma = 0.5

secondLast (x:y:[]) = x
secondLast (_:xs) = secondLast xs

replaceLast xs x = init xs ++ [x]

{-hessianFromSimplex :: (V.Vector Double -> Double) -> [Int] -> [((Int, Int), Double)] -> Simplex -> (V.Vector Double, Matrix Double)
hessianFromSimplex f isInt fixed sim =
  let mat :: [V.Vector Double]
      mat = toRows $ fromColumns $ map fst sim
      fsw ((y0, ymin),ymax) = (y0, max (ymax-y0) (y0-ymin))
      swings = flip map mat $ runStat (fmap fsw $ meanF `both` minFrom 1e80 `both` maxFrom (-1e80)) . toList
      n = length swings
      xv = fromList $ map fst swings
      fxv = f  xv
      fixedpts = map fst fixed
      iswings i | i `elem` isInt = atLeastOne $snd $ swings!!i
                | otherwise = snd $ swings!!i
      funits d i | d/=i = 0
                 | i `elem` isInt = atLeastOne $ snd $ swings!!i
                 | otherwise = snd $ swings!!i
      units = flip map [0..n-1] $ \d -> V.generate n $ funits d
      --http://www.caspur.it/risorse/softappl/doc/sas_docs/ormp/chap5/sect28.htm
      fhess ij@ (i,j) | ij `elem` fixedpts = fromJust $ lookup ij fixed
                      | i>=j =
                         ((f $ xv + units!!i + units!!j)
                          - (f $ xv + units!!i - units!!j)
                          - (f $ xv - units!!i + units!!j)
                          + (f $ xv - units!!i - units!!j) )
                          / (4*(iswings i) * (iswings j))
                      | otherwise = 0.0
      hess1= buildMatrix n n fhess
      hess2 = buildMatrix n n $ \(i,j) ->if i>=j then hess1@@>(i,j)
                                                 else hess1@@>(j,i)

  -- we probably  ought to make pos-definite
  -- http://www.mathworks.com/matlabcentral/newsreader/view_thread/103174
  -- posdefify in R etc
  in (fromList (map (fst) swings), hess2) -}

atLeastOne :: Double -> Double
atLeastOne x | isNaN x || isInfinite x = 1.0
             | x < -1.0 || x > 1.0 = realToFrac $ round x
             | x < 0 = -1.0
             | otherwise = 1.0

genInitial :: (V.Vector Double -> Double) -> [Int] -> (Int -> Double) -> V.Vector Double -> Simplex
genInitial f isInt h x0 = sim where
  n = length $ V.toList x0
  unit d = V.generate n $ \j -> if j /=d then 0.0 else if d `elem` isInt then atLeastOne $ h j*(x0!d)
                                                                          else h j*(x0!d)
  mkv d = with f $ x0 + unit d
  sim = (x0, f x0) : map mkv [0..n-1]
  (!) = (V.!)

goNm :: (V.Vector Double -> Double) -> [Int] -> Double -> Int -> Int -> Simplex -> Simplex
goNm f' isInt tol nmin nmax sim' = go f' 0 $ sortBy (comparing snd) sim'  where
  go f i sim = let nsim = sortBy (comparing snd) $ (nmStep f isInt sim)
                   fdiff = abs $ snd (last nsim) - snd (head nsim)
               in case () of
                      _ |  (fdiff < tol && i>nmin) || i>nmax -> nsim
                        |  all (<0) (map snd sim) && any (>0) (map snd nsim) -> sim
                        |  any (isNaN) (map snd nsim) -> sim
                        |  otherwise   -> go f (i+1) nsim

goNmVerbose :: (V.Vector Double -> Double) -> [Int] -> Double -> Int -> Int -> Simplex -> Simplex
goNmVerbose f' isInt tol nmin nmax sim' = go f' 0 $ sortBy (comparing snd) sim' where
  go f i sim = let nsim = sortBy (comparing snd) $ (nmStep f isInt sim)
                   fdiff = trace ("1: "++ show (fst (head nsim)) ++ "\nlast: "++
                                show (fst (last nsim)) ++ "\n#"++show i++": "++
                                show (map snd nsim))
                            abs $ snd (last nsim) - snd (head nsim)
               in case () of
                    _ |  (fdiff < tol && i>nmin) || i>nmax -> nsim
                      |  all (<0) (map snd sim) && any (>0) (map snd nsim) -> sim
                      |  any (isNaN) (map snd nsim) -> sim
                      |  otherwise   -> go f (i+1) nsim

nmStep :: (V.Vector Double -> Double) -> [Int] -> Simplex -> Simplex
nmStep f isInt s0 = snext where
   x0 = centroid $ init s0
   xnp1 = fst (last s0)
   fxnp1 = snd (last s0)
   xr = x0 + nmAlpha * (x0 - xnp1)
   fxr = f xr
   fx1 = snd $ head s0
   snext = if fx1 <= fxr && fxr <= (snd $ secondLast s0)
              then replaceLast s0 (xr,fxr)
              else sexpand
   xe = x0 + nmGamma * (x0-xnp1)
   fxe = f xe
   sexpand = if fxr > fx1
                then scontract
                else if fxe < fxr
                        then replaceLast s0 (xe,fxe)
                        else replaceLast s0 (xr,fxr)
   xc = xnp1 + scale nmRho  (x0-xnp1)
   fxc = f xc
   scontract = if fxc < fxnp1
                  then replaceLast s0 (xc,fxc)
                  else sreduce
   sreduce = case s0 of
              p0@(x1,_):rest -> p0 : (flip map rest $ \(xi,_) -> with f $ x1+ scale nmRho (xi-x1))

with f x = (x, f x)

max_stars_repo_name: glutamate/probably-base

max_issues_repo_name: glutamate/probably-base

max_forks_repo_name: glutamate/probably-base


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Main (main) where

import Control.Monad.IO.Class (liftIO)
import Data.Complex (Complex)
import Data.Foldable (toList)
import Data.List (nub,
                  sort)
import Text.PrettyPrint.Mainland
import Text.PrettyPrint.Mainland.Class
import Text.Read (readMaybe)

import Spiral
import Spiral.Backend.C
import Spiral.Config
import Spiral.Driver
import Spiral.Exp
import Spiral.Search.OpCount
import Spiral.OpCount
import Spiral.SPL
import Spiral.SPL.Run
import Spiral.Util.Uniq

parseSizesArg :: Monad m => String -> m [Int]
parseSizesArg = fmap concat . mapM parseRange . splitOn (== ',')

parseRange :: Monad m => String -> m [Int]
parseRange s =
    case splitOn (== '-') s of
      [x]   -> do start <- parseInt x
                  return [start]
      [x,y] -> do start <- parseInt x
                  end   <- parseInt y
                  return [start..end]
      _     -> fail $ "Cannot parse range: " ++ s

parseInt :: Monad m => String -> m Int
parseInt s = case readMaybe s of
                 Nothing -> fail $ "Cannot parse integer: " ++ s
                 Just n  -> return n

splitOn :: (Char -> Bool) -> String -> [String]
splitOn p s = case dropWhile p s of
                  "" -> []
                  s' -> w : splitOn p s''
                        where (w, s'') = break p s'

main :: IO ()
main =
    defaultMainWith mempty $ \args -> do
    sizes <- nub . sort . concat <$> mapM parseSizesArg args
    useComplexType <- asksConfig $ testDynFlag UseComplex
    liftIO $ putDocLn $ text "\"size\",\"totalops\",\"mulops\",\"addops\""
    mapM_ (search useComplexType) sizes
  where
    search :: Bool -> Int -> Spiral ()
    search True  n = searchOpCount (DFT n :: SPL (Exp (Complex Double))) >>= go n
    search False n = searchOpCount (Re (DFT n) :: SPL (Exp Double)) >>= go n

    go :: (Typed a, Num (Exp a)) => Int -> SPL (Exp a) -> Spiral ()
    go n e = do
        -- pprint e
        prog <- toProgram ("hspiral_dft_" ++ show n) e
        -- pprint prog
        ops <- countProgramOps prog
        resetUnique
        defs <- evalCg $ cgProgram prog
        outp <- asksConfig output
        case outp of
          Nothing -> return ()
          Just{}  -> writeOutput (toList defs)
        liftIO $ putDocLn $ folddoc (\x y -> x <> comma <> y) [ppr n, ppr (allOps ops), ppr (mulOps ops), ppr (addOps ops)]

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE RoleAnnotations #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ParallelListComp #-}
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE PartialTypeSignatures #-}
{-# OPTIONS_GHC -fno-cse -fno-full-laziness #-}

module Data.Discrimination.Grouping
  ( Group(..)
  , Grouping(..)
  , Grouping1(..)
  -- * Combinators
  , nub, nubWith
  , group, groupWith
  , groupingEq
  , runGroup
  -- * Internals
  , hashing
  , word8s
  ) where

import Control.Monad hiding (mapM_)
import Control.Monad.Primitive
import Control.Monad.ST
import Control.Monad.ST.Unsafe
import Data.Complex
import Data.Discrimination.Internal.WordMap as WordMap
import Data.Foldable hiding (concat)
import Data.Functor.Compose
import Data.Functor.Contravariant
import Data.Functor.Contravariant.Divisible
import Data.Functor.Contravariant.Generic
import Data.Hashable
import Data.Int
import Data.Semigroup hiding (Any)
import Data.Primitive.MutVar
import Data.Primitive.PrimArray
import Data.Promise
import Data.Proxy
import Data.Ratio
import Data.Typeable
import Data.Void
import Data.Word
import GHC.Integer.GMP.Internals
import GHC.Word
import Numeric.Natural
import Prelude hiding (read, concat, mapM_)

-- | Productive Stable Unordered Discriminator

newtype Group a = Group
  { getGroup :: forall m b. PrimMonad m
             => (b -> m (b -> m ())) -> m (a -> b -> m ())
  } deriving Typeable

instance Contravariant Group where
  contramap f m = Group $ \k -> do
    g <- getGroup m k
    return (g . f)

instance Divisible Group where
  conquer = Group $ \ (k :: b -> m (b -> m ())) -> do
    v <- newMutVar undefined
    writeMutVar v $ \b -> k b >>= writeMutVar v
    return $ \ _ b -> readMutVar v >>= ($ b)

  divide f m n = Group $ \k -> do
    kbcd <- getGroup m $ \ (c, d) -> do
      kcd <- getGroup n k
      kcd c d
      return $ uncurry kcd
    return $ \ a d -> case f a of
      (b, c) -> kbcd b (c, d)

instance Decidable Group where
  choose f m n = Group $ \k -> do
    kb <- getGroup m k
    kc <- getGroup n k
    return (either kb kc . f)

  lose k = Group $ \_ -> return (absurd . k)

instance Semigroup (Group a) where
  (<>) = divide (\a -> (a,a))

instance Monoid (Group a) where
  mempty = conquer
  mappend = (<>)

--------------------------------------------------------------------------------
-- Primitives
--------------------------------------------------------------------------------

groupingWord64 :: Group Word64
groupingWord64 = Group $ \k -> do
  mt <- newMutVar WordMap.empty
  return $ \a b -> readMutVar mt >>= \m -> case WordMap.lookup a m of
    Nothing -> k b >>= \p -> writeMutVar mt (insert a p m)
    Just n -> n b

-- | This may be useful for pragmatically accelerating a grouping structure by
-- preclassifying by a hash function
--
-- Semantically,
--
-- @
-- grouping = hashing <> grouping
-- @
hashing :: Hashable a => Group a
hashing = contramap hash grouping

--------------------------------------------------------------------------------
-- * Unordered Discrimination (for partitioning)
--------------------------------------------------------------------------------

-- | 'Eq' equipped with a compatible stable unordered discriminator.
class Grouping a where
  -- | For every surjection @f@,
  --
  -- @
  -- 'contramap' f 'grouping' ≡ 'grouping'
  -- @

  grouping :: Group a
#ifndef HLINT
  default grouping :: Deciding Grouping a => Group a
  grouping = deciding (Proxy :: Proxy Grouping) grouping
#endif

instance Grouping Void where grouping = lose id
instance Grouping () where grouping = conquer
instance Grouping Word8 where grouping = contramap fromIntegral groupingWord64
instance Grouping Word16 where grouping = contramap fromIntegral groupingWord64
instance Grouping Word32 where grouping = contramap fromIntegral groupingWord64
instance Grouping Word64 where grouping = groupingWord64
instance Grouping Word where grouping = contramap fromIntegral groupingWord64
instance Grouping Int8 where grouping = contramap fromIntegral groupingWord64
instance Grouping Int16 where grouping = contramap fromIntegral groupingWord64
instance Grouping Int32 where grouping = contramap fromIntegral groupingWord64
instance Grouping Int64 where grouping = contramap fromIntegral groupingWord64
instance Grouping Int where grouping = contramap fromIntegral groupingWord64
instance Grouping Char where grouping = contramap (fromIntegral . fromEnum) groupingWord64

instance Grouping Bool
instance (Grouping a, Grouping b) => Grouping (a, b)
instance (Grouping a, Grouping b, Grouping c) => Grouping (a, b, c)
instance (Grouping a, Grouping b, Grouping c, Grouping d) => Grouping (a, b, c, d)
instance Grouping a => Grouping [a]
instance Grouping a => Grouping (Maybe a)
instance (Grouping a, Grouping b) => Grouping (Either a b)
instance Grouping a => Grouping (Complex a) where
  grouping = divide (\(a :+ b) -> (a, b)) grouping grouping

instance Grouping Integer where
  grouping = contramap word8s grouping

word8s :: Integer -> [Word8]
word8s i = runST $ unsafeIOToST $ do
  p@(MutablePrimArray mba) :: MutablePrimArray RealWorld Word8 <- newPrimArray (fromIntegral $ W# (sizeInBaseInteger i 256#))
  _ <- exportIntegerToMutableByteArray i mba 0## 1#
  primArrayToList <$> unsafeFreezePrimArray p

instance Grouping Natural where grouping = contramap toInteger grouping

#if __GLASGOW_HASKELL__ >= 800
instance Grouping a => Grouping (Ratio a) where
#else
instance (Grouping a, Integral a) => Grouping (Ratio a) where
#endif
  grouping = divide (\r -> (numerator r, denominator r)) grouping grouping

instance (Grouping1 f, Grouping1 g, Grouping a) => Grouping (Compose f g a) where
  grouping = getCompose `contramap` grouping1 (grouping1 grouping)

class Grouping1 f where
  grouping1 :: Group a -> Group (f a)
#ifndef HLINT
  default grouping1 :: Deciding1 Grouping f => Group a -> Group (f a)
  grouping1 = deciding1 (Proxy :: Proxy Grouping) grouping
#endif

instance Grouping1 []
instance Grouping1 Maybe
instance Grouping a => Grouping1 (Either a)
instance Grouping a => Grouping1 ((,) a)
instance (Grouping a, Grouping b) => Grouping1 ((,,) a b)
instance (Grouping a, Grouping b, Grouping c) => Grouping1 ((,,,) a b c)
instance (Grouping1 f, Grouping1 g) => Grouping1 (Compose f g) where
  grouping1 f = getCompose `contramap` grouping1 (grouping1 f)
instance Grouping1 Complex where
  grouping1 f = divide (\(a :+ b) -> (a, b)) f f

-- | Valid definition for @('==')@ in terms of 'Grouping'.
groupingEq :: Grouping a => a -> a -> Bool
groupingEq a b = runST $ do
  rn <- newMutVar (0 :: Word8)
  k <- getGroup grouping $ \_ -> do
    modifyMutVar' rn (+1)
    return return
  k a ()
  k b ()
  n <- readMutVar rn
  return $ n == 1
{-# INLINE groupingEq #-}

runGroup :: Group a -> [(a,b)] -> [[b]]
runGroup (Group m) xs = runLazy (\p0 -> do
    rp <- newMutVar p0
    f <- m $ \ b -> do
      p <- readMutVar rp
      q <- promise []
      p' <- promise []
      p != (b : demand q) : demand p'
      writeMutVar rp p'
      rq <- newMutVar q
      return $ \b' -> do
        q' <- readMutVar rq
        q'' <- promise []
        q' != b' : demand q''
        writeMutVar rq q''
    mapM_ (uncurry f) xs
  ) []

--------------------------------------------------------------------------------
-- * Combinators
--------------------------------------------------------------------------------

-- | /O(n)/. Similar to 'Data.List.group', except we do not require groups to be clustered.
--
-- This combinator still operates in linear time, at the expense of storing history.
--
-- The result equivalence classes are __not__ sorted, but the grouping is stable.
--
-- @
-- 'group' = 'groupWith' 'id'
-- @
group :: Grouping a => [a] -> [[a]]
group as = runGroup grouping [(a, a) | a <- as]

-- | /O(n)/. This is a replacement for 'GHC.Exts.groupWith' using discrimination.
--
-- The result equivalence classes are __not__ sorted, but the grouping is stable.
groupWith :: Grouping b => (a -> b) -> [a] -> [[a]]
groupWith f as = runGroup grouping [(f a, a) | a <- as]

-- | /O(n)/. This upgrades 'Data.List.nub' from @Data.List@ from /O(n^2)/ to /O(n)/ by using
-- productive unordered discrimination.
--
-- @
-- 'nub' = 'nubWith' 'id'
-- 'nub' as = 'head' 'Control.Applicative.<$>' 'group' as
-- @
nub :: Grouping a => [a] -> [a]
nub = nubWith id

-- | /O(n)/. Online 'nub' with a Schwartzian transform.
--
-- @
-- 'nubWith' f as = 'head' 'Control.Applicative.<$>' 'groupWith' f as
-- @
nubWith :: Grouping b => (a -> b) -> [a] -> [a]
nubWith f xs = runLazy (\p0 -> do
    rp <- newMutVar p0
    k <- getGroup grouping $ \a -> do
      p' <- promise []
      p <- readMutVar rp
      p != a : demand p'
      writeMutVar rp p'
      return $ \ _ -> return ()
    mapM_ (\x -> k (f x) x) xs
  ) []

max_stars_repo_name: alexfmpe/discrimination

max_issues_repo_name: alexfmpe/discrimination

max_forks_repo_name: alexfmpe/discrimination


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}

--
-- Copyright (c) 2009-2011, ERICSSON AB
-- All rights reserved.
--
-- Redistribution and use in source and binary forms, with or without
-- modification, are permitted provided that the following conditions are met:
--
--     * Redistributions of source code must retain the above copyright notice,
--       this list of conditions and the following disclaimer.
--     * Redistributions in binary form must reproduce the above copyright
--       notice, this list of conditions and the following disclaimer in the
--       documentation and/or other materials provided with the distribution.
--     * Neither the name of the ERICSSON AB nor the names of its contributors
--       may be used to endorse or promote products derived from this software
--       without specific prior written permission.
--
-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
-- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-- SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-- CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-- OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-- OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
--

module Feldspar.Core.Constructs.Num
    ( NUM (..)
    ) where

import Language.Syntactic
import Language.Syntactic.Constructs.Binding

import Data.Complex (Complex(..))

import Feldspar.Range
import Feldspar.Core.Types
import Feldspar.Core.Interpretation
import Feldspar.Core.Constructs.Literal
import Feldspar.Core.Constructs.Integral
import Feldspar.Core.Constructs.Complex

data NUM a
  where
    Abs  :: (Type a, Num a, Num (Size a)) => NUM (a :-> Full a)
    Sign :: (Type a, Num a, Num (Size a)) => NUM (a :-> Full a)
    Add  :: (Type a, Num a, Num (Size a)) => NUM (a :-> a :-> Full a)
    Sub  :: (Type a, Num a, Num (Size a)) => NUM (a :-> a :-> Full a)
    Mul  :: (Type a, Num a, Num (Size a)) => NUM (a :-> a :-> Full a)

instance Semantic NUM
  where
    semantics Abs  = Sem "abs" abs
    semantics Sign = Sem "signum" signum
    semantics Add  = Sem "(+)" (+)
    semantics Sub  = Sem "(-)" (-)
    semantics Mul  = Sem "(*)" (*)

instance Equality NUM where equal = equalDefault; exprHash = exprHashDefault
instance Render   NUM where renderArgs = renderArgsDefault
instance ToTree   NUM
instance Eval     NUM where evaluate = evaluateDefault
instance EvalBind NUM where evalBindSym = evalBindSymDefault
instance Sharable NUM

instance AlphaEq dom dom dom env => AlphaEq NUM NUM dom env
  where
    alphaEqSym = alphaEqSymDefault

instance SizeProp (NUM :|| Type)
  where
    sizeProp (C' Abs)  (WrapFull a :* Nil)               = abs (infoSize a)
    sizeProp (C' Sign) (WrapFull a :* Nil)               = signum (infoSize a)
    sizeProp (C' Add)  (WrapFull a :* WrapFull b :* Nil) = infoSize a + infoSize b
    sizeProp (C' Sub)  (WrapFull a :* WrapFull b :* Nil) = infoSize a - infoSize b
    sizeProp (C' Mul)  (WrapFull a :* WrapFull b :* Nil) = infoSize a * infoSize b

instance ( (NUM     :|| Type) :<: dom
         , (Literal :|| Type) :<: dom
         , (INTEGRAL :|| Type) :<: dom
         , (COMPLEX :|| Type) :<: dom
         , OptimizeSuper dom
         )
      => Optimize (NUM :|| Type) dom
  where
    constructFeatOpt _ (C' Abs) (a :* Nil)
        | RangeSet r <- infoRange (getInfo a)
        , isNatural r
        = return a

    constructFeatOpt _ (C' Sign) (a :* Nil)
        | RangeSet ra <- infoRange (getInfo a)
        , 0 `rangeLess` ra
        = return (literalDecor 1)

    constructFeatOpt _ (C' Sign) (a :* Nil)
        | RangeSet ra <- infoRange (getInfo a)
        , ra `rangeLess` 0
        = return (literalDecor (-1))

    constructFeatOpt opts (C' Add) (a :* b :* Nil)
        | Just 0 <- viewLiteral b = return a
        | Just 0 <- viewLiteral a = return b
        | alphaEq a b = constructFeatOpt opts (c' Mul) (a :* literalDecor 2 :* Nil)

    constructFeatOpt opts s@(C' Add) (a :* (op :$ b :$ c) :* Nil)
        | Just al       <- viewLiteral a
        , Just (C' Add) <- prjF op
        , Just cl       <- viewLiteral c
        = constructFeat opts s (b :* literalDecor (al+cl) :* Nil)

    constructFeatOpt opts s@(C' Add) (a :* (op :$ b :$ c) :* Nil)
        | Just al       <- viewLiteral a
        , Just (C' Sub) <- prjF op
        , Just cl       <- viewLiteral c
        = constructFeat opts s (b :* literalDecor (al-cl) :* Nil)

    constructFeatOpt opts s@(C' Add) ((op :$ a :$ b) :* c :* Nil)
        | Just cl       <- viewLiteral c
        , Just (C' Add) <- prjF op
        , Just bl       <- viewLiteral b
        = constructFeat opts s (a :* literalDecor (bl+cl) :* Nil)

    constructFeatOpt opts s@(C' Add) ((op :$ a :$ b) :* c :* Nil)
        | Just cl       <- viewLiteral c
        , Just (C' Sub) <- prjF op
        , Just bl       <- viewLiteral b
        = constructFeat opts s (a :* literalDecor (cl-bl) :* Nil)

    constructFeatOpt opts (C' Add) ((op1 :$ a :$ b) :* (op2 :$ c :$ d) :* Nil)
        | Just (C' Add) <- prjF op1
        , Just (C' Add) <- prjF op2
        , Just bl       <- viewLiteral b
        , Just dl       <- viewLiteral d
        = do
            ac <- constructFeat opts (c' Add) (a :* c :* Nil)
            constructFeat opts (c' Add) (ac :* literalDecor (bl+dl) :* Nil)

    constructFeatOpt opts (C' Add) ((op1 :$ a :$ b) :* (op2 :$ c :$ d) :* Nil)
        | Just (C' Add) <- prjF op1
        , Just (C' Sub) <- prjF op2
        , alphaEq a c
        , alphaEq b d
        = constructFeat opts (c' Add) (a :* c :* Nil)

    -- x `mod` y + y * (x `div` y) ==> x
    -- Partial index calculations materialized from contractT . expandT 2
    -- in MultiDim.hs, which is a no-op.
    constructFeatOpt opts (C' Add) ((rem :$ a :$ b) :* (mul :$ c :$ (quot :$ d :$ e)) :* Nil)
        | Just (C' Rem)  <- prjF rem
        , Just (C' Mul)  <- prjF mul
        , Just (C' Quot) <- prjF quot
        , alphaEq a d
        , alphaEq c e
        , alphaEq b e
        = return a

    -- literal a - (b + literal c) ==> literal (a-c) - b
    -- constructFeatOpt opts s@(C' Sub) (a :* (op :$ b :$ c) :* Nil)
    --     | Just a'      <- viewLiteral a
    --     , Just (C' Add) <- prjF op
    --     , Just c'      <- viewLiteral c
    --     = constructFeat opts s (literalDecor (a'-c') :* b :* Nil)

    -- literal a - (b - literal c) ==> literal (a+c) - b
    -- constructFeatOpt opts s@(C' Sub) (a :* (op :$ b :$ c) :* Nil)
    --     | Just a'      <- viewLiteral a
    --     , Just (C' Sub) <- prjF op
    --     , Just c'      <- viewLiteral c
    --     = constructFeat opts s (literalDecor (a'+c') :* b :* Nil)

    -- (a + literal b) - literal c ==> a + literal (b - c)
    constructFeatOpt opts (C' Sub) ((op :$ a :$ b) :* c :* Nil)
        | Just cl         <- viewLiteral c
        , Just s@(C' Add) <- prjF op
        , Just bl         <- viewLiteral b
        = constructFeat opts s (a :* literalDecor (bl-cl) :* Nil)

    -- (a - literal b) - literal c ==> a - literal (b + c)
    constructFeatOpt opts s@(C' Sub) ((op :$ a :$ b) :* c :* Nil)
        | Just cl       <- viewLiteral c
        , Just (C' Sub) <- prjF op
        , Just bl       <- viewLiteral b
        = constructFeat opts s (a :* literalDecor (bl+cl) :* Nil)

    constructFeatOpt opts (C' Sub) ((op1 :$ a :$ b) :* (op2 :$ c :$ d) :* Nil)
        | Just (C' Add) <- prjF op1
        , Just (C' Sub) <- prjF op2
        , alphaEq a c
        , alphaEq b d
        = constructFeat opts (c' Add) (b :* d :* Nil)

    constructFeatOpt _ (C' Sub) (a :* b :* Nil)
        | Just 0 <- viewLiteral b = return a
        | alphaEq a b             = return $ literalDecor 0

    -- (x + yi) * i ==> -y + xi; (x + yi) * (-i) ==> y - xi
    constructFeatOpt opts (C' Mul) (a :* iunit :* Nil)
        | ComplexType FloatType <- infoType (getInfo iunit)
        , Just (0 :+ k)         <- viewLiteral iunit
        , abs k == 1
        = do
             ra <- constructFeat opts (c' RealPart) (a :* Nil)
             ia <- constructFeat opts (c' ImagPart) (a :* Nil)
             iainv <- constructFeatOpt opts (c' Mul) (literalDecor (-k) :* ia :* Nil)
             rainv <- constructFeatOpt opts (c' Mul) (literalDecor k :* ra :* Nil)
             constructFeatOpt opts (c' MkComplex) (iainv :* rainv :* Nil)

    constructFeatOpt _ (C' Mul) (a :* b :* Nil)
        | Just 0 <- viewLiteral a = return a
        | Just 1 <- viewLiteral a = return b
        | Just 0 <- viewLiteral b = return b
        | Just 1 <- viewLiteral b = return a

    constructFeatOpt opts s@(C' Mul) (a :* (op :$ b :$ c) :* Nil)
        | Just al       <- viewLiteral a
        , Just (C' Mul) <- prjF op
        , Just cl       <- viewLiteral c
        = constructFeat opts s (b :* literalDecor (al*cl) :* Nil)

    constructFeatOpt opts s@(C' Mul) ((op :$ a :$ b) :* c :* Nil)
        | Just cl       <- viewLiteral c
        , Just (C' Mul) <- prjF op
        , Just bl       <- viewLiteral b
        = constructFeat opts s (a :* literalDecor (bl*cl) :* Nil)

    constructFeatOpt opts (C' Mul) ((op1 :$ a :$ b) :* (op2 :$ c :$ d) :* Nil)
        | Just (C' Mul) <- prjF op1
        , Just (C' Mul) <- prjF op2
        , Just b'      <- viewLiteral b
        , Just d'      <- viewLiteral d
        = do
            ac <- constructFeat opts (c' Mul) (a :* c :* Nil)
            constructFeat opts (c' Mul) (ac :* literalDecor (b'*d') :* Nil)

    -- Cases to make sure literals end up to the right:
    constructFeatOpt opts (C' Add) (a :* b :* Nil)
        | Just _ <- viewLiteral a = constructFeatUnOpt opts (c' Add) (b :* a :* Nil)

    constructFeatOpt opts (C' Mul) (a :* b :* Nil)
        | Just _ <- viewLiteral a = constructFeatUnOpt opts (c' Mul) (b :* a :* Nil)

    constructFeatOpt opts a args = constructFeatUnOpt opts a args

    constructFeatUnOpt opts x@(C' _) = constructFeatUnOptDefault opts x

max_stars_repo_name: kffaxen/feldspar-language

max_issues_repo_name: kffaxen/feldspar-language

max_forks_repo_name: kffaxen/feldspar-language


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE MagicHash           #-}
{-# LANGUAGE PatternSynonyms     #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE UnboxedTuples       #-}
{-# LANGUAGE ViewPatterns        #-}
module Numeric.Quaternion.Internal
    ( Quaternion (..), Quater(Quater)
    ) where

import Numeric.DataFrame.Type  (KnownBackend)
import Numeric.Matrix.Internal (Matrix)
import Numeric.PrimBytes       (PrimBytes)
import Numeric.Vector.Internal (Vector)
import Text.Read

-- | @(x,y,z,w)@ of a quaternion
pattern Quater :: Quaternion t => t -> t -> t -> t -> Quater t
pattern Quater a b c d <- (unpackQ# -> (# a, b, c, d #))
  where
    Quater = packQ
{-# COMPLETE Quater #-}

-- | Quaternion operations
class ( Floating (Quater t), Floating t, Ord t, PrimBytes t
      , KnownBackend t '[3], KnownBackend t '[4]
      , KnownBackend t '[3,3], KnownBackend t '[4, 4]
      ) => Quaternion t where
    -- | Quaternion data type. The ordering of coordinates is @(x,y,z,w)@,
    --   where @w@ is the argument, and @x y z@ are the components of a 3D vector
    data Quater t
    -- | Set the quaternion in format @(x,y,z,w)@
    packQ :: t -> t -> t -> t -> Quater t
    -- | Get the values of the quaternion in format @(x,y,z,w)@
    unpackQ# :: Quater t -> (# t, t, t, t #)
    -- | Set the quaternion from 3D axis vector and argument
    fromVecNum :: Vector t 3 -> t -> Quater t
    -- | Set the quaternion from 4D vector in format @(x,y,z,w)@
    fromVec4 :: Vector t 4 -> Quater t
    -- | Transform the quaternion to 4D vector in format @(x,y,z,w)@
    toVec4 :: Quater t -> Vector t 4
    -- | Get scalar square of the quaternion.
    --
    --   >>> realToFrac (square q) == q * conjugate q
    square :: Quater t -> t
    -- | Imaginary part of the quaternion (orientation vector)
    im :: Quater t -> Quater t
    -- | Real part of the quaternion
    re :: Quater t -> Quater t
    -- | Imaginary part of the quaternion as a 3D vector
    imVec :: Quater t -> Vector t 3
    -- | Real part of the quaternion as a scalar
    taker :: Quater t -> t
    -- | i-th component
    takei :: Quater t -> t
    -- | j-th component
    takej :: Quater t -> t
    -- | k-th component
    takek :: Quater t -> t
    -- | Conjugate quaternion (negate imaginary part)
    conjugate :: Quater t -> Quater t
    -- | Rotates and scales vector in 3D using quaternion.
    --   Let \( q = c (\cos \frac{\alpha}{2}, v \sin \frac{\alpha}{2}) \)
    --     , \( c > 0 \), \( {|v|}^2 = 1 \);
    --   then the rotation angle is \( \alpha \), and the axis of rotation is \(v\).
    --   Scaling is proportional to \( c^2 \).
    --
    --   >>> rotScale q x == q * x * (conjugate q)
    rotScale :: Quater t -> Vector t 3 -> Vector t 3
    -- | Creates a quaternion @q@ from two vectors @a@ and @b@,
    --   such that @rotScale q a == b@.
    getRotScale :: Vector t 3 -> Vector t 3 -> Quater t
    -- | Creates a rotation versor from an axis vector and an angle in radians.
    --   Result is always a unit quaternion (versor).
    --   If the argument vector is zero, then result is a real unit quaternion.
    axisRotation :: Vector t 3 -> t -> Quater t
    -- | Quaternion rotation angle \( \alpha \)
    --   (where \( q = c (\cos \frac{\alpha}{2},  v \sin \frac{\alpha}{2}) \)
    --     , \( c > 0 \), \( {|v|}^2 = 1 \)).
    --
    --   >>> q /= 0 ==> axisRotation (imVec q) (qArg q) == signum q
    qArg :: Quater t -> t
    -- | Create a quaternion from a rotation matrix.
    --   Note, that rotations of \(q\) and \(-q\) are equivalent, there result of this
    --   function may be ambiguious. Assume the sign of the result to be chosen arbitrarily.
    fromMatrix33 :: Matrix t 3 3 -> Quater t
    -- | Create a quaternion from a homogenious coordinates trasform matrix.
    --   Ignores matrix translation transform.
    --   Note, that rotations of \(q\) and \(-q\) are equivalent, there result of this
    --   function may be ambiguious. Assume the sign of the result to be chosen arbitrarily.
    fromMatrix44 :: Matrix t 4 4 -> Quater t
    -- | Create a rotation matrix from a quaternion.
    --   Note, that rotations of \(q\) and \(-q\) are equivalent, so the following property holds:
    --
    --   >>> toMatrix33 q == toMatrix33 (-q)
    toMatrix33 :: Quater t -> Matrix t 3 3
    -- | Create a homogenious coordinates trasform matrix from a quaternion.
    --   Translation of the output matrix is zero.
    --   Note, that rotations of \(q\) and \(-q\) are equivalent, so the following property holds:
    --
    --   >>> toMatrix44 q == toMatrix44 (-q)
    toMatrix44 :: Quater t -> Matrix t 4 4

instance (Show t, Quaternion t, Ord t, Num t) => Show (Quater t) where
    showsPrec p (Quater x y z w)
        = case finS of
            SEmpty -> showChar '0'
            Simple -> finF
            SParen -> showParen (p > 6) finF
      where
        (finS, finF) = go SEmpty
          [(w, Nothing), (x, Just 'i'), (y, Just 'j'), (z, Just 'k')]
        go :: ShowState -> [(t, Maybe Char)] -> (ShowState, ShowS)
        go s ((v,l):xs)
          | (s0, f0) <- showComponent s v l
          , (s', f') <- go s0 xs
            = (s', f0 . f')
        go s [] = (s, id)
        showLabel Nothing  = id
        showLabel (Just c) = showChar c
        showComponent :: ShowState -> t -> Maybe Char -> (ShowState, ShowS)
        showComponent sState val mLabel = case (sState, compare val 0) of
          (_     , EQ) -> ( sState, id )
          (SEmpty, GT) -> ( Simple, shows val . showLabel mLabel )
          (SEmpty, LT) -> ( SParen, shows val . showLabel mLabel )
          (_     , GT) -> ( SParen
                          , showString " + " . shows val . showLabel mLabel )
          (_     , LT) -> ( SParen
                          , showString " - " . shows (negate val) . showLabel mLabel )

data ShowState = SEmpty | Simple | SParen
    deriving Eq

instance (Read t, Quaternion t, Num t) => Read (Quater t) where
    readPrec     = parens $ readPrec >>= go id 0 0 0 0
      where
        go :: (t -> t) -> t -> t -> t -> t -> t -> ReadPrec (Quater t)
        go f x y z w new =
          let def = pure (Quater x y z (f new))
              withLabel EOF         = def
              withLabel (Ident "i")
                = (lexP >>= proceed (f new) y z w) <++ pure (Quater (f new) y z w)
              withLabel (Ident "j")
                = (lexP >>= proceed x (f new) z w) <++ pure (Quater x (f new) z w)
              withLabel (Ident "k")
                = (lexP >>= proceed x y (f new) w) <++ pure (Quater x y (f new) w)
              withLabel l           = proceed x y z (f new) l
          in (lexP >>= withLabel) <++ def
        proceed :: t -> t -> t -> t -> Lexeme -> ReadPrec (Quater t)
        proceed x y z w (Symbol "+") = readPrec >>= go id x y z w
        proceed x y z w (Symbol "-") = readPrec >>= go negate x y z w
        proceed x y z w EOF          = pure (Quater x y z w)
        proceed _ _ _ _ _            = pfail

    readListPrec = readListPrecDefault
    readList     = readListDefault

max_stars_repo_name: achirkin/fasttensor

max_issues_repo_name: achirkin/fasttensor

max_forks_repo_name: achirkin/fasttensor


text:
--
-- (c) Susumu Katayama
--
{-# LANGUAGE UndecidableInstances, OverlappingInstances, TemplateHaskell, CPP, FlexibleInstances #-}
-- x #define TESTEQ

module MagicHaskeller.Classification -- (
                      -- randomTestFilter, -- ::  Filtrable a => (b->a) -> Matrix b -> Matrix b
                      -- )
                      where
import Prelude hiding ((/))
#ifdef TFRANDOM
import System.Random.TF.Gen
#else
import System.Random
#endif
import MagicHaskeller.MyCheck
import Data.Char
import Data.List
import Control.Monad
import Control.Monad.Search.Combinatorial
import Data.Complex

import MagicHaskeller.MHTH
import MagicHaskeller.TimeOut
import MagicHaskeller.T10

import MagicHaskeller.Instantiate(compareRealFloat)

class (Search m) => SStrategy m where
    sfilter :: Relation r =>
               (k->k->r) -> (Int->Int) -> m ([k],e) -> m ([k],e)
    ofilter :: Relation r =>
               (k->k->r) -> m (k,e) -> m (k,e)

instance SStrategy Matrix where
    sfilter = sfilterMx
    ofilter = ofilterMx

instance SStrategy DBound where
    sfilter = sfilterDB
    ofilter = ofilterDB

#ifdef TFRANDOM
arbitraries :: Arbitrary a => [a]
arbitraries = arbs 0 (seedTFGen (12279932681387497184,218462391894233257,1625759680792809304,12756118543158863164))
arbs :: Arbitrary a => Int -> TFGen -> [a]
arbs n stdgen  =  case map (splitn stdgen 8) [0..255] of
                    g0:gs -> zipWith f [n..] gs ++ arbs (n+255) g0 -- I think 255 seeds is enough, but just in case.
    where Gen f = arbitrary
#else
arbitraries :: Arbitrary a => [a]
arbitraries = arbs 0 (mkStdGen 1)
arbs :: Arbitrary a => Int -> StdGen -> [a]
arbs n stdgen  =  case split stdgen of
                    (g0,g1) -> f n g0 : arbs (n+1) g1
    where Gen f = arbitrary
#endif

(/~) :: [a] -> (a->a->Bool) -> [[a]]
[]      /~  eq  =  []
(x:xs)  /~  eq  =  case partition (x `eq`) xs of
                     (same, diff) -> (x:same) : (diff /~ eq)
{-
$B$J$*!$>e5-(B(/)$B$*$h$S2<5-(Bnub$B$N>l9g!$L58B%j%9%H$G$b$$$1$k!%(B
*T10> let {[]     / eq = []; (x:xs) / eq = case partition (x `eq`) xs of (same, diff) -> (x:same) : (diff / eq)}
*T10> cycle "hogeha" / (==)
["hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh
(snip)
hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhInterrupted.
*T10> map head $ cycle "hogeha" / (==)
"hogeaInterrupted.
*T10> Data.List.nub $ cycle "hogeha"
"hogeaTerminated
(snip)
*T10> Data.List.nub $ repeat 'h'
"hTerminated
-}

-- ToDo: deal with timeout

-- complete set of representatives
nubSortBy :: (a->a->Ordering) -> [a] -> [a]
nubSortBy = mergesortWithBy const
nubSortByBot :: (a->a->Maybe Ordering) -> [a] -> [a]
nubSortByBot = mergesortWithByBot const
-- quotient set
(/<) :: [a] -> (a->a->Ordering) -> [[a]]
xs /< cmp  =  mergesortWithBy (++)
                (\x y -> head x `cmp` head y)
                (map return xs)
(/<?) :: [a] -> (a->a->Maybe Ordering) -> [[a]]
xs /<? cmp  =  mergesortWithByBot (++)
                (\x y -> head x `cmp` head y)
                (map return xs)
class Eq rel => Relation rel where
 -- remove duplicates, and sort if |rel==Ordering|
    fromListBy   :: (k->k->rel) -> [k] -> [k]
    fromListBy cmp = map head . (/cmp)
-- used to pick up the shallowest expression in |DBound|
    fromListByDB :: (k->k->rel) -> [(k,Int)] -> [(k,Int)]
    fromListByDB rel ts =
          map  (minimumBy (\x y -> compare (snd y) (snd x)))
               (ts / (\x y -> rel (fst x) (fst y)))
    -- \NB : |maximumBy| returns the last of the maxima,
    -- while |minimumBy| the first of the minima.
 -- quotient set
    (/)          :: [k] -> (k->k->rel) -> [[k]]
 -- merge two lists,
    appendWithBy :: (k->k->k) -> -- combiner
                    (k->k->rel) -> [k] -> [k] -> [k]
    diffBy       :: (k -> k -> rel) -> [k] -> [k] -> [k]
 -- counterpart of EQ
    cEQ          :: rel
-- merge two quotient sets
appendQuotientsBy ::
  (Relation rel) =>
  (k -> k -> rel) -> [[k]] -> [[k]] -> [[k]]
appendQuotientsBy rel =
   appendWithBy (++) (\ (x:_) (y:_) -> x `rel` y)
-- merge two complete sets of representatives
appendRepresentativesBy ::
   (Relation rel) =>
   (k -> k -> rel) -> [k] -> [k] -> [k]
appendRepresentativesBy = appendWithBy const

instance Relation Bool where
    fromListBy    = Data.List.nubBy
    (/)           = (/~)
    appendWithBy  = unionWithBy
    diffBy        = Data.List.deleteFirstsBy
    cEQ           = True
unionWithBy combiner eq []      ys = ys
unionWithBy combiner eq (x:xs)  ys =
   case break (eq x) ys of
      (_,   [])    ->  x  :  unionWithBy combiner eq xs ys
      (ts,  h:ds)  ->  (x `combiner` h)
                          :  unionWithBy combiner eq xs (ts++ds)
instance Relation Ordering where
    fromListBy    = nubSortBy
    fromListByDB rel =
        mergesortWithBy
                    (\x y -> if snd x < snd y then y else x)
                    (\x y -> fst x `rel` fst y)
    (/)           = (/<)
    appendWithBy  = mergeWithBy
    diffBy        = diffSortedBy
    cEQ           = EQ
instance Relation (Maybe Ordering) where
    fromListBy    = nubSortByBot
    fromListByDB rel =
        mergesortWithByBot
                    (\x y -> if snd x < snd y then y else x)
                    (\x y -> fst x `rel` fst y)
    (/)           = (/<?)
    appendWithBy  = mergeWithByBot
    diffBy        = diffSortedByBot
    cEQ           = Just EQ

randomTestFilter ::  (SStrategy m, Filtrable a) =>
                     (Int->Int) -> m (e,a) -> m (e,a)
randomTestFilter numRandoms = filt numRandoms . fmap (\ t@(_,a) -> (a,t))

unsafeRandomTestFilter ::  (SStrategy m, Filtrable a) =>
                           Maybe Int -- ^ microsecs until timeout
                               -> (Int->Int) -> m (e,a) -> m (e,a)
unsafeRandomTestFilter mto numRandoms = unsafeFilt mto numRandoms . fmap (\ t@(_,a) -> (a,t))

mapFst f (a,b) = (f a, b)

class Filtrable a where
    filt    :: SStrategy m => (Int->Int) -> m (a,e) -> m e
    filtFun :: (SStrategy m, Arbitrary b) =>
               (Int->Int) -> m (b->a,e) -> m e
    unsafeFilt    :: SStrategy m =>
                     Maybe Int -> (Int->Int) -> m (a,e) -> m e
    unsafeFiltFun :: (SStrategy m, Arbitrary b) =>
                     Maybe Int -> (Int->Int) -> m (b->a,e) -> m e

instance  (Arbitrary a, Filtrable r) => Filtrable (a->r)
  where
    filt      = filtFun
    filtFun f = filt f . fmap (mapFst uncurry)
    unsafeFilt    mto f = unsafeFiltFun mto f
    unsafeFiltFun mto f = unsafeFilt mto f . fmap (mapFst uncurry)

#ifdef TESTEQ
instance Eq a => Filtrable a where
    filt     = filtNullary  (==)
    filtFun  = filtUnary    (==)
#else
instance Ord a => Filtrable a where
    filt     = filtNullary  compare
    filtFun  = filtUnary    compare
    unsafeFilt    mto = filtNullary  (unsafeOpWithPTO mto compare)
    unsafeFiltFun mto = filtUnary    (unsafeOpWithPTO mto compare)
instance Filtrable Double where
    filt     = filtNullary  compareRealFloat
    filtFun  = filtUnary    compareRealFloat
    unsafeFilt    mto = filtNullary  (unsafeOpWithPTO mto compareRealFloat)
    unsafeFiltFun mto = filtUnary    (unsafeOpWithPTO mto compareRealFloat)
#endif
filtNullary ::  (SStrategy m, Relation r) =>
                (k->k->r) -> (Int->Int) -> m (k,e) -> m e
filtNullary  op _ =  fmap snd . ofilter op
filtUnary    op f =  fmap snd . sfilter op f .
                      fmap (mapFst (flip map arbitraries))

instance  (RealFloat a, Ord a) =>
          Filtrable (Complex a) where
    filt     = filtNullary  compareCx
    filtFun  = filtUnary    compareCx
    unsafeFilt    mto = filtNullary  (unsafeOpWithPTO mto compareCx)
    unsafeFiltFun mto = filtUnary    (unsafeOpWithPTO mto compareCx)
compareCx ::  (RealFloat a, Ord a) =>
              Complex a -> Complex a -> Ordering
(a:+b) `compareCx` (c:+d) = case compare a c of
                              EQ -> compare b d
                              o  -> o
ofilterMx ::  Relation r =>
            (k->k->r) -> Matrix (k,e) -> Matrix (k,e)
ofilterMx op (Mx xss)
        = let
            (k,_) `rel` (l,_) = k `op` l
            mapped     = map (fromListBy rel) xss
            cumulative =
                 scanl  (appendRepresentativesBy rel)
                                 [] mapped
          in Mx $ zipWith (diffBy rel) mapped cumulative

ofilterDB :: Relation rel =>
              (k->k->rel) ->
               DBound (k,e) -> DBound (k,e)
ofilterDB cmp (DB f) = DB $
    \n -> fromListByDB  (\(k,_) (l,_) -> cmp k l)
                        (f n)

cumulativeRepresentatives ::
    Relation rel =>
   [a->a->rel] -> Matrix a -> Matrix a
cumulativeRepresentatives relations mx =
    fmap head (cumulativeQuotients relations mx)

representatives ::
    Relation rel =>
   [a->a->rel] -> Matrix a -> Matrix a
representatives relations mx =
   unscanlByList relations $
        cumulativeRepresentatives relations mx
unscanlByList ::  Relation r =>
                  [k->k->r] -> Matrix k -> Matrix k
unscanlByList (_:rels) (Mx (yss@(xs:xss))) =
    Mx $ xs : zipWith3 diffBy rels xss yss

sfilterMx ::  Relation r =>
              (k->k->r) ->
              (Int->Int) ->
              Matrix ([k],e) -> Matrix ([k],e)
sfilterMx rel numRands = representatives (map (liftRelation rel . numRands) [0..])
liftRelation ::  Relation r =>
                 (k->k->r) ->
                    Int -> ([k],e) -> ([k],e) -> r
liftRelation rel len (xs,_) (ys,_) = liftRel rel len xs ys
liftRel _   0   _      _      = cEQ
liftRel _   _   []     []     = cEQ
liftRel rel len (x:xs) (y:ys) =
    case rel x y of
           c  | c == cEQ   -> liftRel rel (len-1) xs ys
              | otherwise  -> c

sfilterDB ::  Relation rel =>
               (k->k->rel) ->
              (Int->Int)->
                DBound ([k],e) -> DBound ([k],e)
sfilterDB rel numRands (DB f) = DB $ \n -> fromListByDB  (liftRelation rel (numRands n)) (f n)

cumulativeQuotients relations (Mx xss)
   =  let yss:ysss = zipWith (/) xss relations
      in Mx $ scanl  (\rec (r,z) ->
                       appendQuotientsBy r (rec>>=(/r)) z)
                     yss  (zip (tail relations) ysss)

ns = [6..]

max_stars_repo_name: liconoc/MagicWrangler

max_issues_repo_name: junjihashimoto/MagicHaskeller

max_forks_repo_name: junjihashimoto/MagicHaskeller