packages feed

hanalyze-0.2.0.0: src/Hanalyze/Design/Orthogonal.hs

{-# LANGUAGE OverloadedStrings #-}
-- |
-- Module      : Hanalyze.Design.Orthogonal
-- Description : 直交表 (Taguchi 流 @Lₙ@ 表) の標準表・因子割付・出力レンダリング
-- Copyright   : (c) 2026 Aelysce Project (Toshiaki Honda)
-- License     : BSD-3-Clause
--
-- Orthogonal arrays (Taguchi-style @Lₙ@ tables).
--
--   * 'OA'              — orthogonal-array representation (name / run
--     count / factor count / column levels / body).
--   * 'standardArrays'  — standard tables L4 / L8 / L9 / L12 / L16 / L18.
--   * 'lookupOA'        — fetch a standard table by name (e.g. @\"L9\"@).
--   * 'assignFactors'   — bind factors and level values.
--   * 'renderCSV' / 'renderTSV' / 'renderPretty' — emit the run table.
--
-- Two-level series (L8, L16, ...) are generated by @mkL2k@. L4 / L9 /
-- L12 / L18 are defined manually (Plackett-Burman and mixed-level
-- arrays are not derivable from simple subset products).
module Hanalyze.Design.Orthogonal
  ( -- * 型
    OA (..)
  , LevelValue (..)
  , FactorSpec (..)
  , AssignedDesign (..)
    -- * Standard arrays
  , l4
  , l8
  , l9
  , l12
  , l16
  , l18
  , l27
  , standardArrays
  , lookupOA
  , listArrays
  , OAMetadata (..)
  , listArraysWithSize
    -- * 2-level array generation
  , mkL2k
    -- * Factor assignment
  , assignFactors
    -- * Rendering
  , renderRawCSV
  , renderRawTSV
  , renderRawPretty
  , renderCSV
  , renderTSV
  , renderPretty
  ) where

import Data.Bits (testBit, popCount, (.&.), bit)
import Data.Text (Text)
import qualified Data.Text as T
import Text.Printf (printf)

-- ---------------------------------------------------------------------------
-- 型
-- ---------------------------------------------------------------------------

-- | An orthogonal array. Stored as a @runs × cols@ table of 1-based
-- level codes.
data OA = OA
  { oaName    :: Text     -- ^ Display name, e.g. @\"L9(3^4)\"@.
  , oaRuns    :: Int      -- ^ Number of runs.
  , oaFactors :: Int      -- ^ Maximum number of factors (= columns).
  , oaLevels  :: [Int]    -- ^ Level count per column (length 'oaFactors').
  , oaTable   :: [[Int]]  -- ^ Body of the table (@runs × cols@) of
                          --   1-based level codes.
  } deriving (Show, Eq)

-- | A factor level value (text or numeric).
data LevelValue = LText Text | LNumeric Double
  deriving (Show, Eq)

-- | User-supplied factor: a name plus a list of level values.
data FactorSpec = FactorSpec
  { fsName   :: Text
  , fsLevels :: [LevelValue]
  } deriving (Show, Eq)

-- | A run table after factor assignment.
data AssignedDesign = AssignedDesign
  { adArray   :: OA
  , adFactors :: [FactorSpec]
  , adRows    :: [[LevelValue]]
  } deriving (Show, Eq)

-- ---------------------------------------------------------------------------
-- 標準表 (手動定義)
-- ---------------------------------------------------------------------------

-- | L4(2³) — 4 runs, up to 3 two-level factors.
l4 :: OA
l4 = OA "L4(2^3)" 4 3 (replicate 3 2)
  [ [1,1,1]
  , [1,2,2]
  , [2,1,2]
  , [2,2,1]
  ]

-- | L9(3⁴) — 9 runs, up to 4 three-level factors.
l9 :: OA
l9 = OA "L9(3^4)" 9 4 (replicate 4 3)
  [ [1,1,1,1]
  , [1,2,2,2]
  , [1,3,3,3]
  , [2,1,2,3]
  , [2,2,3,1]
  , [2,3,1,2]
  , [3,1,3,2]
  , [3,2,1,3]
  , [3,3,2,1]
  ]

-- | L12(2¹¹) — 12 runs, up to 11 two-level factors (Plackett-Burman).
-- Main effects only (interactions are distributed across all columns).
l12 :: OA
l12 = OA "L12(2^11)" 12 11 (replicate 11 2)
  [ [1,1,1,1,1,1,1,1,1,1,1]
  , [1,1,1,1,1,2,2,2,2,2,2]
  , [1,1,2,2,2,1,1,1,2,2,2]
  , [1,2,1,2,2,1,2,2,1,1,2]
  , [1,2,2,1,2,2,1,2,1,2,1]
  , [1,2,2,2,1,2,2,1,2,1,1]
  , [2,1,2,2,1,1,2,2,1,2,1]
  , [2,1,2,1,2,2,2,1,1,1,2]
  , [2,1,1,2,2,2,1,2,2,1,1]
  , [2,2,2,1,1,1,1,2,2,1,2]
  , [2,2,1,2,1,2,1,1,1,2,2]
  , [2,2,1,1,2,1,2,1,2,2,1]
  ]

-- | L18(2¹×3⁷) — 18 runs, up to 8 factors (column 1 has 2 levels, the
-- remaining 7 columns each have 3 levels).
--
-- One of the most recommended Taguchi-style arrays; can measure main
-- effects plus the column-1 × column-2 interaction.
l18 :: OA
l18 = OA "L18(2^1*3^7)" 18 8 (2 : replicate 7 3)
  [ [1,1,1,1,1,1,1,1]
  , [1,1,2,2,2,2,2,2]
  , [1,1,3,3,3,3,3,3]
  , [1,2,1,1,2,2,3,3]
  , [1,2,2,2,3,3,1,1]
  , [1,2,3,3,1,1,2,2]
  , [1,3,1,2,1,3,2,3]
  , [1,3,2,3,2,1,3,1]
  , [1,3,3,1,3,2,1,2]
  , [2,1,1,3,3,2,2,1]
  , [2,1,2,1,1,3,3,2]
  , [2,1,3,2,2,1,1,3]
  , [2,2,1,2,3,1,3,2]
  , [2,2,2,3,1,2,1,3]
  , [2,2,3,1,2,3,2,1]
  , [2,3,1,3,2,3,1,2]
  , [2,3,2,1,3,1,2,3]
  , [2,3,3,2,1,2,3,1]
  ]

-- | L27(3¹³) — 27 runs, up to 13 three-level factors.
--
-- GF(3) 上の線形形式で生成する (手写しの転記ミスを避ける)。 27 run を
-- @(a,b,c) ∈ {0,1,2}³@ で添字し、 13 列は 3 次元 GF(3) の相異なる 1 次元部分空間
-- を代表する係数ベクトル (先頭非零 = 1 に正規化) で @α·a+β·b+γ·c mod 3@ を取る。
-- 相異なる 1 次元部分空間ゆえ任意の 2 列は強度 2 直交 (各水準組が均等出現)。
-- 列順は標準タグチ L27 の配置 (col1,2 = 基本、 3,4 = その交互作用、 5 = 第3基本…)。
l27 :: OA
l27 = OA "L27(3^13)" 27 13 (replicate 13 3) table
  where
    coeffs :: [(Int, Int, Int)]
    coeffs =
      [ (1,0,0), (0,1,0), (1,1,0), (1,2,0), (0,0,1), (1,0,1), (1,0,2)
      , (0,1,1), (0,1,2), (1,1,1), (1,1,2), (1,2,1), (1,2,2) ]
    table =
      [ [ 1 + ((al*a + be*b + ga*c) `mod` 3) | (al, be, ga) <- coeffs ]
      | r <- [0 .. 26 :: Int]
      , let a = r `div` 9
            b = (r `div` 3) `mod` 3
            c = r `mod` 3 ]

-- ---------------------------------------------------------------------------
-- 2 水準系の生成
-- ---------------------------------------------------------------------------

-- | Build @L_{2^k}(2^{2^k − 1})@ in Taguchi's standard column ordering
-- (column @j@'s value is
-- popCount(j ∧ revBits k r) のパリティ)。
mkL2k :: Int -> OA
mkL2k k =
  OA
    { oaName    = T.pack ("L" ++ show n ++ "(2^" ++ show m ++ ")")
    , oaRuns    = n
    , oaFactors = m
    , oaLevels  = replicate m 2
    , oaTable   = [ [ levelAt r j | j <- [1 .. m] ] | r <- [0 .. n - 1] ]
    }
  where
    n = 2 ^ k
    m = n - 1
    -- Taguchi の標準的な列ラベル順 (col 1 は最上位ビット相当) に合わせるため
    -- 行インデックスをビット反転する。
    revBits :: Int -> Int
    revBits r = sum [ if testBit r i then bit (k - 1 - i) else 0
                    | i <- [0 .. k - 1] ]
    levelAt r j = 1 + (popCount (j .&. revBits r) `mod` 2)

-- | L8(2⁷) — 8 runs, up to 7 two-level factors (generated).
l8 :: OA
l8 = mkL2k 3

-- | L16(2¹⁵) — 16 runs, up to 15 two-level factors (generated).
l16 :: OA
l16 = mkL2k 4

-- ---------------------------------------------------------------------------
-- ルックアップ
-- ---------------------------------------------------------------------------

-- | The standard arrays bundled with the library.
standardArrays :: [OA]
standardArrays = [l4, l8, l9, l12, l16, l18, l27]

-- | Look up a standard array by short name (e.g. @\"L9\"@).
lookupOA :: Text -> Maybe OA
lookupOA name0 = case T.toUpper name0 of
  "L4"  -> Just l4
  "L8"  -> Just l8
  "L9"  -> Just l9
  "L12" -> Just l12
  "L16" -> Just l16
  "L18" -> Just l18
  "L27" -> Just l27
  _     -> Nothing

-- | List of available orthogonal arrays (used by CLI @doe list@).
listArrays :: [(Text, Text)]
listArrays = [ (oaName a, descr a) | a <- standardArrays ]
  where
    descr a =
      T.pack (show (oaRuns a)) <> " runs, max "
      <> T.pack (show (oaFactors a)) <> " factors"

-- | Structured metadata for an orthogonal array. Suitable for UI
-- listings that want to filter / sort by run count or level pattern.
data OAMetadata = OAMetadata
  { omName    :: !Text   -- ^ e.g. @\"L9(3^4)\"@.
  , omRuns    :: !Int    -- ^ Number of runs.
  , omFactors :: !Int    -- ^ Maximum number of factors.
  , omLevels  :: ![Int]  -- ^ Level count per column.
  , omDescr   :: !Text   -- ^ Free-form description (matches 'listArrays').
  } deriving (Show, Eq)

-- | Same coverage as 'listArrays' but with structured fields.
listArraysWithSize :: [OAMetadata]
listArraysWithSize =
  [ OAMetadata (oaName a) (oaRuns a) (oaFactors a) (oaLevels a)
               (T.pack (show (oaRuns a)) <> " runs, max "
                <> T.pack (show (oaFactors a)) <> " factors")
  | a <- standardArrays
  ]

-- ---------------------------------------------------------------------------
-- 因子割当
-- ---------------------------------------------------------------------------

-- | Assign user-supplied factor names and level values to the columns of
-- an orthogonal array, returning the expanded run table.
--
-- - 因子数が表の列数を超えるとエラー
-- - 各因子の水準数が割当先列の水準数と一致しないとエラー
assignFactors :: OA -> [FactorSpec] -> Either Text AssignedDesign
assignFactors oa specs
  | nSpecs > oaFactors oa =
      Left $ "Too many factors: " <> oaName oa
             <> " has only " <> T.pack (show (oaFactors oa)) <> " columns; got "
             <> T.pack (show nSpecs)
  | not (null mismatches) =
      Left $ "Factor level mismatch: " <> T.intercalate "; " mismatches
  | otherwise =
      Right AssignedDesign
        { adArray   = oa
        , adFactors = specs
        , adRows    = [ [ fsLevels (specs !! (j - 1)) !! (lvl - 1)
                        | (j, lvl) <- zip [1 .. nSpecs] (take nSpecs row) ]
                      | row <- oaTable oa ]
        }
  where
    nSpecs    = length specs
    expected  = take nSpecs (oaLevels oa)
    actuals   = map (length . fsLevels) specs
    mismatches =
      [ fsName (specs !! i) <> " expected " <> T.pack (show e)
        <> " levels, got " <> T.pack (show a)
      | (i, (e, a)) <- zip [0..] (zip expected actuals)
      , e /= a ]

-- ---------------------------------------------------------------------------
-- 出力
-- ---------------------------------------------------------------------------

-- | Render an orthogonal array as raw CSV (columns are @F1, F2, …@).
renderRawCSV :: OA -> Text
renderRawCSV oa = renderRawWith "," oa

-- | Render an orthogonal array as raw TSV.
renderRawTSV :: OA -> Text
renderRawTSV oa = renderRawWith "\t" oa

-- | Render an orthogonal array as a delimiter-separated table.
renderRawWith :: Text -> OA -> Text
renderRawWith sep oa =
  let header = T.intercalate sep
                 [ "F" <> T.pack (show j) | j <- [1 .. oaFactors oa] ]
      body   = T.intercalate "\n"
                 [ T.intercalate sep [ T.pack (show v) | v <- row ]
                 | row <- oaTable oa ]
  in header <> "\n" <> body <> "\n"

-- | Pretty-print a named orthogonal array with aligned columns.
renderRawPretty :: OA -> Text
renderRawPretty oa =
  let names    = "Run" : [ "F" <> T.pack (show j) | j <- [1 .. oaFactors oa] ]
      colWidth = maximum (map T.length names) `max` 3
      pad t    = let n = colWidth - T.length t
                 in T.replicate n " " <> t
      header   = T.intercalate "  " (map pad names)
      body     = T.intercalate "\n"
                   [ T.intercalate "  "
                       (pad (T.pack (show r))
                       : [ pad (T.pack (show v)) | v <- row ])
                   | (r, row) <- zip [1::Int ..] (oaTable oa) ]
  in T.pack (T.unpack (oaName oa)) <> "\n" <> header <> "\n" <> body

-- | Render a factor-assigned run table as CSV.
renderCSV :: AssignedDesign -> Text
renderCSV = renderWith ","

-- | Render a factor-assigned run table as TSV.
renderTSV :: AssignedDesign -> Text
renderTSV = renderWith "\t"

-- | Render a factor-assigned run table with a custom field separator.
renderWith :: Text -> AssignedDesign -> Text
renderWith sep ad =
  let header = T.intercalate sep ("Run" : map fsName (adFactors ad))
      body   = T.intercalate "\n"
                 [ T.intercalate sep (T.pack (show r) : map fmtLevel row)
                 | (r, row) <- zip [1::Int ..] (adRows ad) ]
  in header <> "\n" <> body <> "\n"

fmtLevel :: LevelValue -> Text
fmtLevel (LText t)    = t
fmtLevel (LNumeric d)
  | d == fromIntegral (round d :: Integer) = T.pack (show (round d :: Integer))
  | otherwise                              = T.pack (printf "%g" d)

-- | Pretty-print a factor-assigned run table.
renderPretty :: AssignedDesign -> Text
renderPretty ad =
  let names      = "Run" : map fsName (adFactors ad)
      cells      =
        [ T.pack (show r) : map fmtLevel row
        | (r, row) <- zip [1::Int ..] (adRows ad) ]
      colWidths  = map (\i -> maximum (map (T.length . safeIx i)
                                       (names : cells)))
                       [0 .. length names - 1]
      safeIx i xs = if i < length xs then xs !! i else ""
      pad i t    = let n = colWidths !! i - T.length t
                   in T.replicate n " " <> t
      fmtRow row = T.intercalate "  "
                     [ pad i (safeIx i row) | i <- [0 .. length names - 1] ]
  in oaName (adArray ad) <> "  (" <> T.pack (show (oaRuns (adArray ad)))
     <> " runs, " <> T.pack (show (length (adFactors ad)))
     <> " of " <> T.pack (show (oaFactors (adArray ad))) <> " columns assigned)\n"
     <> fmtRow names <> "\n"
     <> T.intercalate "\n" (map fmtRow cells)