rubic

module AdjAverage where

import Debug.Trace
import KdTree
import Data.List
import Data.Maybe
import ReturnNearestMap
import qualified Data.Array.Repa as R
import Basis

rubicCube =
        let
           
p0 = Point3dIndexed (Point3d -2 -2 -2) 0 [(Point3d -2 -2 -3),(Point3d -2 -2 -1),(Point3d -2 -3 -2),(Point3d -1 -2 -2),(Point3d -2 -1 -2),(Point3d -3 -2 -2)]
p1 = Point3dIndexed (Point3d -2 -2 0) 1 [(Point3d -2 -2 -1),(Point3d -2 -2 1),(Point3d -2 -3 0),(Point3d -1 -2 0),(Point3d -2 -1 0),(Point3d -3 -2 0)]
p2 = Point3dIndexed (Point3d -2 -2 2) 2 [(Point3d -2 -2 1),(Point3d -2 -2 3),(Point3d -2 -3 2),(Point3d -1 -2 2),(Point3d -2 -1 2),(Point3d -3 -2 2)]
p3 = Point3dIndexed (Point3d -2 0 -2) 3 [(Point3d -2 0 -3),(Point3d -2 0 -1),(Point3d -2 -1 -2),(Point3d -1 0 -2),(Point3d -2 1 -2),(Point3d -3 0 -2)]
p4 = Point3dIndexed (Point3d -2 0 0) 4 [(Point3d -2 0 -1),(Point3d -2 0 1),(Point3d -2 -1 0),(Point3d -1 0 0),(Point3d -2 1 0),(Point3d -3 0 0)]
p5 = Point3dIndexed (Point3d -2 0 2) 5 [(Point3d -2 0 1),(Point3d -2 0 3),(Point3d -2 -1 2),(Point3d -1 0 2),(Point3d -2 1 2),(Point3d -3 0 2)]
p6 = Point3dIndexed (Point3d -2 2 -2) 6 [(Point3d -2 2 -3),(Point3d -2 2 -1),(Point3d -2 1 -2),(Point3d -1 2 -2),(Point3d -2 3 -2),(Point3d -3 2 -2)]
p7 = Point3dIndexed (Point3d -2 2 0) 7 [(Point3d -2 2 -1),(Point3d -2 2 1),(Point3d -2 1 0),(Point3d -1 2 0),(Point3d -2 3 0),(Point3d -3 2 0)]
p8 = Point3dIndexed (Point3d -2 2 2) 8 [(Point3d -2 2 1),(Point3d -2 2 3),(Point3d -2 1 2),(Point3d -1 2 2),(Point3d -2 3 2),(Point3d -3 2 2)]
p9 = Point3dIndexed (Point3d 0 -2 -2) 9 [(Point3d 0 -2 -3),(Point3d 0 -2 -1),(Point3d 0 -3 -2),(Point3d 1 -2 -2),(Point3d 0 -1 -2),(Point3d -1 -2 -2)]
p10 = Point3dIndexed (Point3d 0 -2 0) 10 [(Point3d 0 -2 -1),(Point3d 0 -2 1),(Point3d 0 -3 0),(Point3d 1 -2 0),(Point3d 0 -1 0),(Point3d -1 -2 0)]
p11 = Point3dIndexed (Point3d 0 -2 2) 11 [(Point3d 0 -2 1),(Point3d 0 -2 3),(Point3d 0 -3 2),(Point3d 1 -2 2),(Point3d 0 -1 2),(Point3d -1 -2 2)]
p12 = Point3dIndexed (Point3d 0 0 -2) 12 [(Point3d 0 0 -3),(Point3d 0 0 -1),(Point3d 0 -1 -2),(Point3d 1 0 -2),(Point3d 0 1 -2),(Point3d -1 0 -2)]
p13 = Point3dIndexed (Point3d 0 0 0) 13 [(Point3d 0 0 -1),(Point3d 0 0 1),(Point3d 0 -1 0),(Point3d 1 0 0),(Point3d 0 1 0),(Point3d -1 0 0)]
p14 = Point3dIndexed (Point3d 0 0 2) 14 [(Point3d 0 0 1),(Point3d 0 0 3),(Point3d 0 -1 2),(Point3d 1 0 2),(Point3d 0 1 2),(Point3d -1 0 2)]
p15 = Point3dIndexed (Point3d 0 2 -2) 15 [(Point3d 0 2 -3),(Point3d 0 2 -1),(Point3d 0 1 -2),(Point3d 1 2 -2),(Point3d 0 3 -2),(Point3d -1 2 -2)]
p16 = Point3dIndexed (Point3d 0 2 0) 16 [(Point3d 0 2 -1),(Point3d 0 2 1),(Point3d 0 1 0),(Point3d 1 2 0),(Point3d 0 3 0),(Point3d -1 2 0)]
p17 = Point3dIndexed (Point3d 0 2 2) 17 [(Point3d 0 2 1),(Point3d 0 2 3),(Point3d 0 1 2),(Point3d 1 2 2),(Point3d 0 3 2),(Point3d -1 2 2)]
p18 = Point3dIndexed (Point3d 2 -2 -2) 18 [(Point3d 2 -2 -3),(Point3d 2 -2 -1),(Point3d 2 -3 -2),(Point3d 3 -2 -2),(Point3d 2 -1 -2),(Point3d 1 -2 -2)]
p19 = Point3dIndexed (Point3d 2 -2 0) 19 [(Point3d 2 -2 -1),(Point3d 2 -2 1),(Point3d 2 -3 0),(Point3d 3 -2 0),(Point3d 2 -1 0),(Point3d 1 -2 0)]
p20 = Point3dIndexed (Point3d 2 -2 2) 20 [(Point3d 2 -2 1),(Point3d 2 -2 3),(Point3d 2 -3 2),(Point3d 3 -2 2),(Point3d 2 -1 2),(Point3d 1 -2 2)]
p21 = Point3dIndexed (Point3d 2 0 -2) 21 [(Point3d 2 0 -3),(Point3d 2 0 -1),(Point3d 2 -1 -2),(Point3d 3 0 -2),(Point3d 2 1 -2),(Point3d 1 0 -2)]
p22 = Point3dIndexed (Point3d 2 0 0) 22 [(Point3d 2 0 -1),(Point3d 2 0 1),(Point3d 2 -1 0),(Point3d 3 0 0),(Point3d 2 1 0),(Point3d 1 0 0)]
p23 = Point3dIndexed (Point3d 2 0 2) 23 [(Point3d 2 0 1),(Point3d 2 0 3),(Point3d 2 -1 2),(Point3d 3 0 2),(Point3d 2 1 2),(Point3d 1 0 2)]
p24 = Point3dIndexed (Point3d 2 2 -2) 24 [(Point3d 2 2 -3),(Point3d 2 2 -1),(Point3d 2 1 -2),(Point3d 3 2 -2),(Point3d 2 3 -2),(Point3d 1 2 -2)]
p25 = Point3dIndexed (Point3d 2 2 0) 25 [(Point3d 2 2 -1),(Point3d 2 2 1),(Point3d 2 1 0),(Point3d 3 2 0),(Point3d 2 3 0),(Point3d 1 2 0)]
p26 = Point3dIndexed (Point3d 2 2 2) 26 [(Point3d 2 2 1),(Point3d 2 2 3),(Point3d 2 1 2),(Point3d 3 2 2),(Point3d 2 3 2),(Point3d 1 2 2)]


        in
                [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15,p16,p17,p18,p19,p20,p21,p22,p23,p24,p25,p26]


getRelativeTwnPos relv =
    case relv of
        1->(-1,-1,-1)
        2->(1,-1,-1)
        3->(1,1,-1)
        4->(-1,1,-1)
        5->(-1,-1,1)
        6->(1,-1,1)
        7->(1,1,1)
        8->(-1,1,1)
        9->(0,-1,-1)
        10->(1,0,-1)
        11->(0,1,-1)
        12->(-1,0,-1)
        13->(0,-1,1)
        14->(1,0,1)
        15->(0,1,1)
        16->(-1,0,1)
        17->(-1,-1,0)
        18->(1,-1,0)
        19->(1,1,0)
        20->(-1,1,0)

--relv:: 左下から1_indexed の　1~8頂点
-- 返り値::[[(int,int,int,int)]] (x,y,z,1~8頂点) 
getRelVertsNearRelVert relv =
     case relv of
         1 -> [(-1,-1,-1,7),(0,-1,-1,8),(-1,0,-1,6),(0,0,-1,5),(-1,-1,0,3),(0,-1,0,4),(-1,0,0,2),(0,0,0,1)]
         2 -> [(0,-1,-1,7),(1,-1,-1,8),(0,0,-1,6),(1,0,-1,5),(0,-1,0,3),(1,-1,0,4),(0,0,0,2),(1,0,0,1)]
         3 -> [(0,0,-1,7),(1,0,-1,8),(0,1,-1,6),(1,1,-1,5),(0,0,0,3),(1,0,0,4),(0,1,0,2),(1,1,0,1)]
         4 -> [(-1,0,-1,7),(0,0,-1,8),(-1,1,-1,6),(0,1,-1,5),(-1,0,0,3),(0,0,0,4),(-1,1,0,2),(0,1,0,1)]
         5 -> [(-1,-1,0,7),(0,-1,0,8),(-1,0,0,6),(0,0,0,5),(-1,-1,1,3),(0,-1,1,4),(-1,0,1,2),(0,0,1,1)]
         6 -> [(0,-1,0,7),(1,-1,0,8),(0,0,0,6),(1,0,0,5),(0,-1,1,3),(1,-1,1,4),(0,0,1,2),(1,0,1,1)]
         7 -> [(0,0,0,7),(1,0,0,8),(0,1,0,6),(1,1,0,5),(0,0,1,3),(1,0,1,4),(0,1,1,2),(1,1,1,1)]
         8 -> [(-1,0,0,7),(0,0,0,8),(-1,1,0,6),(0,1,0,5),(-1,0,1,3),(0,0,1,4),(-1,1,1,2),(0,1,1,1)]


getCellIndexAndRelVerts:: Int  -> [[Int]] -> [[(Int,Int)]]
getCellIndexAndRelVerts cellindex adjacent_surfs =  
     
     let
         getAdjCellIndex :: (Int,Int,Int) -> Int ->Maybe Int
         getAdjCellIndex (dx,dy,dz) cellindex =
            if cellindex < 0 then Nothing 
            else if dx /= 0 then
                 let num = (if dx == 1 then 3 else 5) in 
                     getAdjCellIndex  (0,dy,dz) ((adjacent_surfs!! cellindex) !! num)
            else if dy /= 0 then
                 let num = (if dy == 1 then 4 else 2) in 
                     getAdjCellIndex  (0,0,dz)  ((adjacent_surfs !! cellindex) !! num)
            else if dz /= 0 then
                 let num = (if dz == 1 then 1 else 0) in 
                     getAdjCellIndex  (0,0,0)  ((adjacent_surfs !! cellindex) !! num)
            else 
                 Just cellindex
     in 
         do
             i <- [1..8]
             let nearverts = getRelVertsNearRelVert i 
             return $ catMaybes $ map (\(dx,dy,dz,relv ) -> case getAdjCellIndex (dx,dy,dz) cellindex of {Just x -> Just (x,relv);Nothing -> Nothing}) nearverts

getCellIndexfromPQR (p,q,r) pmax qmax rmax = (p*qmax*rmax) + (q*rmax) + r

getPQRfromCellIndex index pmax qmax rmax =
     let
          r = index `mod` rmax
          q = ((index - r) `div` rmax) `mod` qmax
          p = (index - (q*rmax) -r) `div` (qmax*rmax)
     in
         (p,q,r)

getAdjVertsfromCellIndex cellindex  pmax qmax rmax =
    let
         (p,q,r) = getPQRfromCellIndex cellindex pmax qmax rmax
    in
         do
             i <- [1..8]
             let vs =  do
                  (dx,dy,dz,v) <- getRelVertsNearRelVert i
                  if 0 <= (p+dx) &&(p+dx) <pmax && 0 <= (q+dy) && (q+dy) < qmax && 0 <= (r+dz) && (r+dz) < rmax then
                       return  $ Just  ((p+dx), (q+dy), (r+dz), v)
                  else
                       return Nothing
             return $ catMaybes vs


getAdjVertsfromPQR (p,q,r) pmax qmax rmax =
    do
        i <- [1..8]
        let vs =  do
             (dx,dy,dz,v) <- getRelVertsNearRelVert i
             if 0 <= (p+dx) &&(p+dx) <pmax && 0 <= (q+dy) && (q+dy) < qmax && 0 <= (r+dz) && (r+dz) < rmax then
                  return  $ Just  (getCellIndexfromPQR (p+dx,q+dy,r+dz) pmax qmax rmax, v)
             else
                  return Nothing
        return $ catMaybes vs

--セル番号から位置関係を特定できる場合に使う
getAdjAverage_old:: (Int,Int,Int) -> Coord R.U  -> Coord R.D
getAdjAverage_old (numr,numt,numz) row  = 
    R.fromFunction (R.Z R.:. numr R.:. numt R.:. numz R.:. 20 R.:. 3) ave
    where
        average xs = realToFrac (sum xs) / fromIntegral (length xs)
        ave v@(R.Z R.:. p R.:. q R.:. r R.:. relv R.:. axis) =  
             let adjs = getAdjVertsfromPQR (p,q,r) numr numt numz in
                 case relv of
                 _ | relv <= 7 -> trace (show (p,q,r,relv) ++  show (adjs !! relv)) $  average (map (\(cellindex,v_) ->(
                                                                    let (p_,q_,r_) = (getPQRfromCellIndex cellindex  numr numt numz) 
                                                                        val        =  row  R.! (R.Z R.:. p_ R.:. q_ R.:. r_ R.:. v_-1 R.:. axis) in trace(show val)val )) (adjs !! relv)) 
                   | otherwise -> row R.! v

--共有頂点の情報を用意した場合に使う
--pos_func:: セルインデックス -> 1_indexed頂点 -> 自然座標(x,y,z)
--relv:: 0~19頂点 
--shared_vertexs (shared_vertexs !! cell) !! num -> num方向の面に隣接するセルの番号 numの並びは(0:-z,1:+z,2:-y,3:+x,4:+y,5:-x)
-- eg: (shared_vetexs !! 10) !! 1 =　10番目のセルの+dz方向の面に隣接するセル番号
-- returnNearestMapの出力をいれる。

getAdjAverage:: Int->Int->(Int->Int -> R.Array R.U R.DIM1 Double) ->[[Int]] -> (Double,Double,Double)
getAdjAverage cellindex relv  pos_func adjacent_surfs   = (ave (R.Z R.:.0),ave (R.Z R.:.1) , ave (R.Z R.:.2))
        where
        average xs = realToFrac (sum xs) / fromIntegral (length xs)
        ave (R.Z R.:. axis) =  
             let adjs = getCellIndexAndRelVerts cellindex adjacent_surfs in
                 case relv of
                 _ | relv <= 7 -> average (map (\(index_,v_) ->let val = (pos_func index_ v_) R.! (R.Z R.:. axis) in  val) (adjs !! relv)) 
                   | otherwise ->trace ("getAdjAverage: relv must be <=7") undefined

--for test 
makeCubeTwn:: R.Array  R.D R.DIM4 Int
makeCubeTwn = R.fromFunction(R.Z R.:.3 R.:. 3 R.:.3 R.:.20) (\( R.Z R.:.x R.:.y R.:.z R.:.n) ->(x*3*3*20)+(y*3*20)+(z*20)+n)

--for test
makeCubeRow:: Coord R.U
makeCubeRow = R.computeS $ R.fromFunction(R.Z R.:.3 R.:. 3 R.:.3 R.:.20 R.:.3) f where
    f (R.Z R.:.p R.:.q R.:.r R.:.n R.:.axis) = 
        let v@(Point3dIndexed (Point3d x y z) _ _) = rubicCube !! (p*9+3*q+r)
            (dx,dy,dz) =getRelativeTwnPos (n+1) in
        case axis of
          0 -> x + dx
          1 -> y + dy
          2 -> z + dz

calcAverage:: Coord R.U ->[[Int]]-> Coord R.U
calcAverage row adjacent_surfs = 
    let 
        sh@(R.Z R.:. numx R.:. numy R.:. numz R.:. numv R.:. numaxis) = R.extent row

        ave (R.Z R.:. x R.:. y R.:. z R.:. relv) =
             if relv <= 7 then
                   getAdjAverage (getCellIndexfromPQR (x,y,z) numx numy numz) (relv) (\x -> \y ->(let (p,q,r) = getPQRfromCellIndex x numx numy numz in R.computeS (R.slice row  (R.Z R.:. p R.:. q R.:.r R.:.y-1 R.:. R.All)))) adjacent_surfs
             else let pos =(( R.computeS (R.slice row (R.Z R.:.x R.:.y R.:.z R.:. relv R.:.R.All)))::R.Array R.U R.DIM1 Double)  in (pos R.! (R.Z R.:.0),pos R.! (R.Z R.:.1),pos R.! (R.Z R.:.2))

        tmp:: R.Array R.U R.DIM4 (Double,Double,Double)
        tmp = R.computeS $ R.fromFunction (R.Z R.:.numx R.:.numy R.:.numz R.:.numv) ave 
    in
       R.computeS $  R.fromFunction sh (\(R.Z R.:.x R.:.y R.:.z R.:.relv R.:.axis) -> (let (px,py,pz) =  tmp R.! (R.Z R.:.x R.:.y R.:.z R.:.relv) in case axis of{ 0-> px; 1->py; 2->pz} )) 


move_p3 :: Point3d->Point3d-> Point3d
move_p3 p d = Point3d (p3x p + p3x d) (p3y p + p3y d) (p3z p + p3z d)

move_p3i :: Point3dIndexed->Point3d->Point3dIndexed
move_p3i p3 d = Point3dIndexed (move_p3 (point p3) d) (index p3) (map (\x -> move_p3 x d) ( surfaceArray p3))

someFunc :: IO ()
someFunc = do
    putStrLn "someFunc"
    --putStrLn $ show $ returnNearestMap [rubicCube !! 3, rubicCube !! 4]
    --let ps' = map (\x -> move_p3i x (Point3d 100 100 100)) rubicCube
    --let ps03 = [ps' !! 3, ps'!! 4]
    --putStrLn $ show $ ps03
    --putStrLn $ show $ returnNearestMap ps03
    --putStrLn $ show  $ (1,2,3)
    --putStrLn $ show  $ getAdjVertsfromPQR (2,0,0) 3 3 3
    print "cube"
    let cube =  makeCubeRow
    let avedcube_old =  ((R.computeS $  getAdjAverage_old (3,3,3) makeCubeRow ) :: Coord R.U)
    print $ cube
    print $ ((R.computeS $(cube R.-^ avedcube_old)) :: Coord R.U)
    --print $ map  (\i -> getCellIndexAndRelVerts i (returnNearestMap rubicCube)) [0..26]
    --print $ makeCubeRow R.! (R.Z R.:. )
    print $ getAdjAverage 9 6 (\x -> \y ->(let (p,q,r) = getPQRfromCellIndex x 3 3 3 in R.computeS (R.slice makeCubeRow (R.Z R.:. p R.:. q R.:.r R.:.y-1 R.:. R.All)))) (returnNearestMap rubicCube)
    let avedcube = calcAverage makeCubeRow (returnNearestMap rubicCube)
    print $ ((R.computeS $ (cube R.-^ avedcube))::Coord R.U)