{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example of GWLP calculation for mixed-level designs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains example code to compute the type-specific generalized word-length pattern (GWLP) for mixed-level designs, more specifically for regular four-and-two-level designs.\n",
    "\n",
    "In four-and-two-level designs, the four-level factors are constructed us the grouping scheme of [Wu & Zhang (1989)](https://doi.org/10.1093/biomet/80.1.203) where the levels of a four-level factor $A$ are based on the levels of three two-level factors $a_1$, $a_2$ and $a_3$, called the pseudo-factors, where $I=a_1a_2a_3$:\n",
    "$$\n",
    "\\begin{array}{ccccc}\n",
    "        a_1 & a_2 & a_3 & & A \\\\\n",
    "        1 & 1 & 0 & \\rightarrow & 0 \\\\\n",
    "        1 & 0 & 1 & \\rightarrow & 1 \\\\\n",
    "        0 & 1 & 1 & \\rightarrow & 2 \\\\\n",
    "        0 & 0 & 0 & \\rightarrow & 3 \\\\\n",
    "    \\end{array}\n",
    "$$\n",
    "The two-level factors can either be main factors or be entirely aliased with a combination of the main factors. Such aliased factors are called added factors. If one of the main factors, used in an added factor, is also used as pseudo-factor in a four-level factor, then the added factor has type $I$. If it is used as pseudo-factor in two distinct four-level factor, it has type $II$, etc $\\ldots$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import oapackage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a $4^{1}2^{7}$ four-and-two-level regular design in 32 runs with 1 four-level factor and 7 two-level factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exampleArray 56: design in OA(32, 42^{18})\n",
      "0000000\n",
      "0100110\n",
      "0101010\n",
      "0001100\n",
      "0110011\n",
      "0010101\n",
      "0011001\n",
      "0111111\n",
      "1010100\n",
      "1110010\n",
      "1111110\n",
      "1011000\n",
      "1100111\n",
      "1000001\n",
      "1001101\n",
      "1101011\n",
      "2101101\n",
      "2001011\n",
      "2000111\n",
      "2100001\n",
      "2011110\n",
      "2111000\n",
      "2110100\n",
      "2010010\n",
      "3111001\n",
      "3011111\n",
      "3010011\n",
      "3110101\n",
      "3001010\n",
      "3101100\n",
      "3100000\n",
      "3000110\n"
     ]
    }
   ],
   "source": [
    "array = oapackage.exampleArray(56, 1)\n",
    "array = array.selectFirstColumns(7)\n",
    "arrayclass = oapackage.arraylink2arraydata(array)\n",
    "array.showarraycompact()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The GWLP for mixed-level designs is computed using the adapted forms of equations (7) and (8) in [Xu and Wu (2001)](https://www.doi.org/10.1214/aos/1009210552).\n",
    "\n",
    "## Distance distribution\n",
    "First, the distance distribution is first computed for all combinations of $(i,j)$ with $i \\in \\{0,1\\}$ and $j=0,\\ldots,7$ and $i\\neq j$. In this example with a $4^12^{7}$ design, that is a $2 \\times 8$ matrix $B$, where $B_{i,j}$ is the number of rows that have $(1-i)$ different elements if their four-level parts and $(7-j)$ different elements in their two-level parts, divided by the number of runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "distance distribution for mixed-level design\n",
      "[[1 0 1 4 1 0 1]\n",
      " [0 2 6 8 6 2 0]]\n"
     ]
    }
   ],
   "source": [
    "Dm = oapackage.distance_distribution_mixed(array, 0)\n",
    "D = np.array(Dm).astype(int)\n",
    "print(f\"distance distribution for mixed-level design\\n{D}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is verified since the sum of the matrix is 32 which is equal to $\\binom{32}{2} 32^{-1}$.\n",
    "\n",
    "## McWilliams Transforms\n",
    "\n",
    "Now, the second step in the computation of the GWLP, is to use the McWilliams Transforms to obtain the GWLP from the distance distribution matrix. The McWilliams Transform will create another $2 \\times 8$ matrix $B^{\\prime}$, obtained using the following formula\n",
    "$$\n",
    "B_{j_{1}, j_{2}}^{\\prime}=N^{-1} \\sum_{i_{1}=0}^{m} \\sum_{i_{2}=0}^{n} B_{i_{1}, i_{2}} P_{j_{1}}\\left(i_{1} ; 1, 4\\right) P_{j_{2}}\\left(i_{2} ; 7, 2\\right)\n",
    "$$\n",
    "where $P_j(x,n,s)$ is the Krawtchouck polynomials for a total of $n$ factors with $s$ levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MacWilliams transform:\n",
      "[[1. 0. 0. 0. 1. 0. 0.]\n",
      " [0. 0. 2. 0. 4. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "N = array.n_rows\n",
    "factor_levels_for_groups = arrayclass.factor_levels_column_groups()\n",
    "Bprime = oapackage.macwilliams_transform_mixed(Dm, N, factor_levels_for_groups, verbose=0)\n",
    "print(\"MacWilliams transform:\")\n",
    "print(np.array(Bprime).astype(int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This matrix is equivalent to the type specific generalized word-length pattern, where the row index indicates the type of words (0 to 1) and the column index indicates the length of the words (0 to 7).\n",
    "\n",
    "The generalized word-length pattern ($A$) can be obtained by summing the rows of the $B^{\\prime}$ matrix anti-diagonally. That is:\n",
    "$$A_j(D) = \\sum_{i^{\\prime}+j^{\\prime}=j} B^{\\prime}_{i^{\\prime},j^{\\prime}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GWLP: (1.0, 0.0, 0.0, 2.0, 1.0, 4.0, 0.0, 0.0)\n"
     ]
    }
   ],
   "source": [
    "gwlp_mixed = oapackage.GWLPmixed(array, 0)\n",
    "print(f\"GWLP: {gwlp_mixed}\")"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "3a893c4079a4533c9db639b850e04910700d284bba4afe56b82dceecf522f906"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}