Skip to content

Instantly share code, notes, and snippets.

@Gitmoko

Gitmoko/rubic

Last active December 30, 2017 12:18
Show Gist options
  • Save Gitmoko/8903fb2a2b789b28fc88a4ce6149e88c to your computer and use it in GitHub Desktop.
Save Gitmoko/8903fb2a2b789b28fc88a4ce6149e88c to your computer and use it in GitHub Desktop.
Download ZIP
Raw
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment