baryonPC_xfang-wfirstxso.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "list_SD36 = ['6x2pt_WFIRSTwidexSO5','6x2pt_WFIRSTstdxSO5']\n",
    "list_COV36 = ['cov_wfirstwidexso_6x2pt.npy','cov_wfirststdxso_6x2pt.npy']\n",
    "\n",
    "list_outfiles = ['6x2pt_WFIRSTwidexSO5','6x2pt_WFIRSTstdxSO5','3x2pt_WFIRSTwidexSO5','3x2pt_WFIRSTstdxSO5','kk_WFIRSTwidexSO5','kk_WFIRSTstdxSO5',\n",
    "                'k3x2pt_WFIRSTwidexSO5','k3x2pt_WFIRSTstdxSO5']\n",
    "\n",
    "dir_datav = \"datav/\"\n",
    "dir_bary = \"datav/\"\n",
    "dir_COV = \"cov/\"\n",
    "dir_LPC = \"baryons/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class gen_LPCi():\n",
    "    def __init__(self,id_SD,outfile,id_include=None):\n",
    "        '''\n",
    "            id_SD = 0~1, int\n",
    "        '''\n",
    "        self.id_SD = id_SD\n",
    "        self.id_include = id_include\n",
    "        self.outfile = outfile\n",
    "        self.list_bary_scenario = ['_eagle','_illustris','_TNG100','_mb2','_owls_AGN','_HzAGN']\n",
    "        self.datav_dmo = self.load_datav(dir_datav,bary_scenario=\"_dmo\",isCutted=1)\n",
    "        print(self.id_SD)\n",
    "        print(self.id_include)\n",
    "        \n",
    "        self.load_COV()\n",
    "        print(self.id_SD)\n",
    "        print(self.datav_dmo.shape, self.COV_cut.shape)\n",
    "        \n",
    "        self.cal_invL()\n",
    "        self.Ndata = len(self.datav_dmo)\n",
    "        self.Nscenario = len(self.list_bary_scenario)\n",
    "        \n",
    "        self.build_Ratio()\n",
    "        self.build_Delta()\n",
    "        self.build_DeltaChy()\n",
    "        self.SVD(in_Delta=self.DeltaChy)\n",
    "        self.write_LPC()\n",
    "        \n",
    "    \n",
    "    def load_datav(self,dir_path,bary_scenario,isCutted=0):\n",
    "        fname_datav = dir_path+list_SD36[self.id_SD]+bary_scenario\n",
    "        print(fname_datav)\n",
    "        datav = np.genfromtxt(fname_datav,skip_header=0,dtype=\"int,double\",usecols=[0,1],names=[\"ind\",\"obs\"])[\"obs\"]\n",
    "        if self.id_include is None:\n",
    "            self.id_include = np.arange(datav.size)\n",
    "        datav = datav[self.id_include]\n",
    "        \n",
    "        self.Ndata_full  = len(datav)\n",
    "        self.takeout_ID  = np.where(datav > 1e-20)[0]\n",
    "        self.zero_ID     = np.where(datav < 1e-20)[0]\n",
    "        \n",
    "        if isCutted==1:\n",
    "            return np.array(datav)[self.takeout_ID]\n",
    "        else:\n",
    "            return np.array(datav)\n",
    "        \n",
    "            \n",
    "    def load_COV(self):\n",
    "        \n",
    "        self.fname_COV = dir_COV+list_COV36[self.id_SD]\n",
    "        print(self.fname_COV)\n",
    "        self.COV_full = np.load(self.fname_COV)\n",
    "        self.COV_full = self.COV_full[self.id_include][:,self.id_include]\n",
    "        \n",
    "        self.COV_cut  = self.COV_full[self.takeout_ID][:,self.takeout_ID]\n",
    "    \n",
    "    \n",
    "    def cal_invL(self):\n",
    "        self.L    = np.linalg.cholesky(self.COV_cut)\n",
    "        self.invL =  np.linalg.inv(self.L)\n",
    "        \n",
    "    def build_Ratio(self):\n",
    "        self.Ratio     = np.zeros((self.Ndata,self.Nscenario))\n",
    "        \n",
    "        for j in range(self.Nscenario):\n",
    "        \n",
    "            datav_bary = self.load_datav(dir_bary,bary_scenario=self.list_bary_scenario[j],isCutted=1)\n",
    "            self.Ratio.T[j] = datav_bary/self.datav_dmo\n",
    "            \n",
    "    def build_Delta(self):\n",
    "        DeltaT = self.Ratio.T*self.datav_dmo-self.datav_dmo\n",
    "        self.Delta = DeltaT.T\n",
    "    \n",
    "    def build_DeltaChy(self):\n",
    "        self.DeltaChy = np.dot(self.invL,self.Delta)\n",
    "        \n",
    "    def SVD(self,in_Delta):\n",
    "        # self.U stores PC modes\n",
    "        # PC1 = self.U.T[0]  ; PC2 = self.U.T[1]\n",
    "        self.U, self.Sdig, VT = np.linalg.svd(in_Delta,full_matrices=True)\n",
    "    \n",
    "    def gen_datav_cosmolike_format(self,datav_cut):\n",
    "        \n",
    "        datav_full = np.zeros(self.Ndata_full)\n",
    "        datav_full[self.takeout_ID] = datav_cut\n",
    "        \n",
    "        return datav_full\n",
    "\n",
    "    \n",
    "    def write_LPC(self):\n",
    "        \n",
    "        self.fname_out = dir_LPC+'LPC_'+self.outfile\n",
    "        \n",
    "        self.LPC1 = self.gen_datav_cosmolike_format(np.dot(self.L,self.U.T[0]))\n",
    "        self.LPC2 = self.gen_datav_cosmolike_format(np.dot(self.L,self.U.T[1]))\n",
    "        self.LPC3 = self.gen_datav_cosmolike_format(np.dot(self.L,self.U.T[2]))\n",
    "        self.LPC4 = self.gen_datav_cosmolike_format(np.dot(self.L,self.U.T[3]))\n",
    "        self.LPC5 = self.gen_datav_cosmolike_format(np.dot(self.L,self.U.T[4]))\n",
    "        self.LPC6 = self.gen_datav_cosmolike_format(np.dot(self.L,self.U.T[5]))\n",
    "        \n",
    "        col_names='LPC1 LPC2 LPC3 LPC4 LPC5 LPC6'\n",
    "        \n",
    "        LPC_Matrix=np.column_stack((self.LPC1,self.LPC2,self.LPC3,self.LPC4,self.LPC5,self.LPC6))\n",
    "\n",
    "        np.savetxt(self.fname_out,LPC_Matrix, fmt='%.8e %.8e %.8e %.8e %.8e %.8e')\n",
    "\n",
    "       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "datav/6x2pt_WFIRSTwidexSO5_dmo\n",
      "0\n",
      "[   0    1    2 ... 1767 1768 1769]\n",
      "cov/cov_wfirstwidexso_6x2pt.npy\n",
      "0\n",
      "(1518,) (1518, 1518)\n",
      "datav/6x2pt_WFIRSTwidexSO5_eagle\n",
      "datav/6x2pt_WFIRSTwidexSO5_illustris\n",
      "datav/6x2pt_WFIRSTwidexSO5_TNG100\n",
      "datav/6x2pt_WFIRSTwidexSO5_mb2\n",
      "datav/6x2pt_WFIRSTwidexSO5_owls_AGN\n",
      "datav/6x2pt_WFIRSTwidexSO5_HzAGN\n"
     ]
    }
   ],
   "source": [
    "PCA = gen_LPCi(0, list_outfiles[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cal_Qexp(scenario):\n",
    "    datav_bary = PCA.load_datav(dir_path=dir_datav, bary_scenario='_'+scenario, isCutted=1)\n",
    "    diff = datav_bary - PCA.datav_dmo\n",
    "    diff_chy = np.dot(PCA.invL, diff)\n",
    "    Qexp = np.dot(PCA.U.T,diff_chy)\n",
    "    return Qexp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "datav/6x2pt_WFIRSTwidexSO5_dmo\n",
      "datav/6x2pt_WFIRSTwidexSO5_eagle\n",
      "datav/6x2pt_WFIRSTwidexSO5_illustris\n",
      "datav/6x2pt_WFIRSTwidexSO5_TNG100\n",
      "datav/6x2pt_WFIRSTwidexSO5_mb2\n",
      "datav/6x2pt_WFIRSTwidexSO5_owls_AGN\n",
      "datav/6x2pt_WFIRSTwidexSO5_HzAGN\n"
     ]
    }
   ],
   "source": [
    "baryons = ['dmo','eagle','illustris','TNG100','mb2','owls_AGN','HzAGN']\n",
    "Qexps = [None]*7\n",
    "for i in range(7):\n",
    "    Qexps[i] = cal_Qexp(scenario=baryons[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dmo [0. 0. 0. 0.]\n",
      "eagle [ 9.3651934   4.18908947  1.17206957 -0.89395842]\n",
      "illustris [75.16754881 -4.71244298  1.58944541 -0.15149137]\n",
      "TNG100 [10.87121777  3.19203493  1.61409652 -0.49579595]\n",
      "mb2 [-1.12450778 -1.99199349 -2.84070904 -0.96812738]\n",
      "owls_AGN [59.92378485  3.41357589 -2.55891895  0.2782546 ]\n",
      "HzAGN [17.68769198  4.15517458  0.12139181  0.4176067 ]\n"
     ]
    }
   ],
   "source": [
    "for i in range(7):\n",
    "    print(baryons[i], (Qexps[i])[0:4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'a' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-15-c2371523eb52>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDeltaChy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'a' is not defined"
     ]
    }
   ],
   "source": [
    "print(a.DeltaChy.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.000e+00 1.000e+00 2.000e+00 ... 1.767e+03 1.768e+03 1.769e+03]\n"
     ]
    }
   ],
   "source": [
    "ells = np.loadtxt(dir_datav+list_SD36[0]+'_dmo')[:,0]\n",
    "print(ells)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: MacOSX\n"
     ]
    }
   ],
   "source": [
    "%matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(6):\n",
    "    plt.plot(np.arange(1518),a.DeltaChy[:,i],label='%d'%(i))\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "datav/6x2pt_WFIRSTwidexSO5_dmo\n",
      "0\n",
      "[   0    1    2 ... 1767 1768 1769]\n",
      "cov/cov_wfirstwidexso_6x2pt.npy\n",
      "0\n",
      "(1518,) (1518, 1518)\n",
      "datav/6x2pt_WFIRSTwidexSO5_eagle\n",
      "datav/6x2pt_WFIRSTwidexSO5_illustris\n",
      "datav/6x2pt_WFIRSTwidexSO5_TNG100\n",
      "datav/6x2pt_WFIRSTwidexSO5_mb2\n",
      "datav/6x2pt_WFIRSTwidexSO5_owls_AGN\n",
      "datav/6x2pt_WFIRSTwidexSO5_HzAGN\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<__main__.gen_LPCi at 0x11080e9d0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gen_LPCi(0, list_outfiles[0])\n",
    "# gen_LPCi(1, list_outfiles[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "datav/6x2pt_WFIRSTwidexSO5_dmo\n",
      "0\n",
      "[   0    1    2 ... 1137 1138 1139]\n",
      "cov/cov_wfirstwidexso_6x2pt.npy\n",
      "0\n",
      "(1012,) (1012, 1012)\n",
      "datav/6x2pt_WFIRSTwidexSO5_eagle\n",
      "datav/6x2pt_WFIRSTwidexSO5_illustris\n",
      "datav/6x2pt_WFIRSTwidexSO5_TNG100\n",
      "datav/6x2pt_WFIRSTwidexSO5_mb2\n",
      "datav/6x2pt_WFIRSTwidexSO5_owls_AGN\n",
      "datav/6x2pt_WFIRSTwidexSO5_HzAGN\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<__main__.gen_LPCi at 0x10c2cb490>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nsource = 10\n",
    "nshear = nsource*(nsource+1)/2\n",
    "ntheta = 15\n",
    "nlens = 10\n",
    "\n",
    "ndata3x2pt= 1455 - (nlens + nsource + 1)*ntheta\n",
    "gen_LPCi(0, list_outfiles[2], np.arange(ndata3x2pt))\n",
    "\n",
    "# ndata3x2pt= 1770 - (nlens + nsource + 1)*ntheta\n",
    "# gen_LPCi(1, list_outfiles[3], np.arange(ndata3x2pt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "datav/6x2pt_WFIRSTwidexSO5_dmo\n",
      "0\n",
      "[1755 1756 1757 1758 1759 1760 1761 1762 1763 1764 1765 1766 1767 1768\n",
      " 1769]\n",
      "cov/cov_wfirstwidexso_6x2pt.npy\n",
      "0\n",
      "(15,) (15, 15)\n",
      "datav/6x2pt_WFIRSTwidexSO5_eagle\n",
      "datav/6x2pt_WFIRSTwidexSO5_illustris\n",
      "datav/6x2pt_WFIRSTwidexSO5_TNG100\n",
      "datav/6x2pt_WFIRSTwidexSO5_mb2\n",
      "datav/6x2pt_WFIRSTwidexSO5_owls_AGN\n",
      "datav/6x2pt_WFIRSTwidexSO5_HzAGN\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<__main__.gen_LPCi at 0x112b71730>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ntheta = 15\n",
    "gen_LPCi(0, list_outfiles[4], 1755+np.arange(ntheta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "datav/6x2pt_WFIRSTwidexSO5_dmo\n",
      "0\n",
      "[1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468\n",
      " 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482\n",
      " 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496\n",
      " 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510\n",
      " 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524\n",
      " 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538\n",
      " 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552\n",
      " 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566\n",
      " 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580\n",
      " 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594\n",
      " 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608\n",
      " 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622\n",
      " 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636\n",
      " 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650\n",
      " 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664\n",
      " 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678\n",
      " 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692\n",
      " 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706\n",
      " 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720\n",
      " 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 1733 1734\n",
      " 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748\n",
      " 1749 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761 1762\n",
      " 1763 1764 1765 1766 1767 1768 1769]\n",
      "cov/cov_wfirstwidexso_6x2pt.npy\n",
      "0\n",
      "(274,) (274, 274)\n",
      "datav/6x2pt_WFIRSTwidexSO5_eagle\n",
      "datav/6x2pt_WFIRSTwidexSO5_illustris\n",
      "datav/6x2pt_WFIRSTwidexSO5_TNG100\n",
      "datav/6x2pt_WFIRSTwidexSO5_mb2\n",
      "datav/6x2pt_WFIRSTwidexSO5_owls_AGN\n",
      "datav/6x2pt_WFIRSTwidexSO5_HzAGN\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<__main__.gen_LPCi at 0x112bc4a90>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ntheta = 15\n",
    "gen_LPCi(0, list_outfiles[6], 1455+np.arange(315)) # k3x2pt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}