062718_R_locfit_19p5.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "locfit 1.5-9.1 \t 2013-03-22\n"
     ]
    }
   ],
   "source": [
    "library(locfit)\n",
    "library(RcppCNPy)\n",
    "library(reticulate)\n",
    "np <- import('numpy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p0 = c(1.11553, 13.06008, 0.44578, 11.29134, 11.75068, 0., 0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##ngal(1)+wp(19)+ggl(20)+vpf(20)+Pcic(40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "seedlist <- c(161,448,712,831,866)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data <- c(np$load('Run_062718/bolp_19p5_um0_b_161.npz'),\n",
    "         np$load('Run_062718/bolp_19p5_um0_b_448.npz'),\n",
    "         np$load('Run_062718/bolp_19p5_um0_b_712.npz'),\n",
    "         np$load('Run_062718/bolp_19p5_um0_b_831.npz'),\n",
    "         np$load('Run_062718/bolp_19p5_um0_b_866.npz'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "setwd('~/Assembly Bias Project/Fisher analysis/Run_062718/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GCV for alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "alphas <- seq(0.1,0.9,0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "calc_gcv_score <- function(i, j){\n",
    "    x <- data[[i]]$f[['param']][(500*j-499):(500*j),j]\n",
    "    gcv_score <- matrix(0,ncol=17,nrow=100)\n",
    "    for (k in seq(100)){\n",
    "        y <- data[[i]]$f[['func_all']][(500*j-499):(500*j),k]\n",
    "        for (l in seq(17)){\n",
    "            gcv_score[k,l] <- as.numeric(gcv(y~x,kern='tria',alpha=alphas[l]))[4]\n",
    "        }\n",
    "    }\n",
    "    gcv_score\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for (i in seq(5)){\n",
    "    np$savez(sprintf('gcv_bolp_20p0_um0_%i.npz',seedlist[i]),gcv_alpha=calc_gcv_score(i,1),\n",
    "            gcv_logM1=calc_gcv_score(i,2),gcv_sigmalogM=calc_gcv_score(i,3),gcv_logM0=calc_gcv_score(i,4),\n",
    "            gcv_logMmin=calc_gcv_score(i,5),gcv_Acen=calc_gcv_score(i,6),gcv_Asat=calc_gcv_score(i,7))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### locfit with alpha=0.3,0.5,0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "setwd('~/Assembly Bias Project/Fisher analysis/Run_062718/bolp_20p0_um0_locfit_test_alpha')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "find_dfdp_test <- function(data, num, alpha10, plot=0){\n",
    "    dfdp = matrix(0,nrow=7,ncol=100)\n",
    "    for (i in seq(7)){\n",
    "        x <- data$f[['param']][(500*i-499):(500*i),i]\n",
    "        p <- p0[i]\n",
    "        for (j in seq(100)){\n",
    "            y <- data$f[['func_all']][(500*i-499):(500*i),j]\n",
    "            fit <- locfit(y~x,deriv=1,alpha=alpha10/10.,kern='tria')\n",
    "            fit_orig <- locfit(y~x,alpha=alpha10/10.,kern='tria')\n",
    "            f <- as.numeric(preplot(fit_orig,p)[2])\n",
    "            fp <- as.numeric(preplot(fit,p)[2])\n",
    "            dfdp[i,j] <- fp\n",
    "            if (plot==1){\n",
    "                png(filename=sprintf('run%i_p%if%i_a%i.png',num,(i-1),(j-1),alpha10))\n",
    "                plot(x, y)\n",
    "                lines(fit_orig,col='blue')\n",
    "                abline(v=p,lty=2)\n",
    "                abline(a=-p*fp+f,b=fp,col='red')\n",
    "                dev.off()\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "    dfdp\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dfdp_236_3 = find_dfdp_test(data[[1]],seedlist[1],alpha10=3,plot=1)\n",
    "#dfdp_483_3 = find_dfdp_test(data[[2]],seedlist[2],alpha10=3,plot=1)\n",
    "#dfdp_518_3 = find_dfdp_test(data[[3]],seedlist[3],alpha10=3,plot=1)\n",
    "#dfdp_563_3 = find_dfdp_test(data[[4]],seedlist[4],alpha10=3,plot=1)\n",
    "#dfdp_831_3 = find_dfdp_test(data[[5]],seedlist[5],alpha10=3,plot=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dfdp_236_5 = find_dfdp_test(data[[1]],seedlist[1],alpha10=5,plot=1)\n",
    "#dfdp_483_5 = find_dfdp_test(data[[2]],seedlist[2],alpha10=5,plot=1)\n",
    "#dfdp_518_5 = find_dfdp_test(data[[3]],seedlist[3],alpha10=5,plot=1)\n",
    "#dfdp_563_5 = find_dfdp_test(data[[4]],seedlist[4],alpha10=5,plot=1)\n",
    "#dfdp_831_5 = find_dfdp_test(data[[5]],seedlist[5],alpha10=5,plot=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dfdp_236_7 = find_dfdp_test(data[[1]],seedlist[1],alpha10=7,plot=1)\n",
    "#dfdp_483_7 = find_dfdp_test(data[[2]],seedlist[2],alpha10=7,plot=1)\n",
    "#dfdp_518_7 = find_dfdp_test(data[[3]],seedlist[3],alpha10=7,plot=1)\n",
    "#dfdp_563_7 = find_dfdp_test(data[[4]],seedlist[4],alpha10=7,plot=1)\n",
    "#dfdp_831_7 = find_dfdp_test(data[[5]],seedlist[5],alpha10=7,plot=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Use 0.7(default) for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "setwd('~/Assembly Bias Project/Fisher analysis/Run_062718/bolp_19p5_um0_locfit_tcub_alpha07/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "find_dfdp <- function(data, num, plot=0){\n",
    "    dfdp = matrix(0,nrow=7,ncol=100)\n",
    "    for (i in seq(7)){\n",
    "        x <- data$f[['param']][(500*i-499):(500*i),i]\n",
    "        p <- p0[i]\n",
    "        for (j in seq(100)){\n",
    "            y <- data$f[['func_all']][(500*i-499):(500*i),j]\n",
    "            fit <- locfit(y~x,deriv=1)\n",
    "            fit_orig <- locfit(y~x)\n",
    "            f <- as.numeric(preplot(fit_orig,p)[2])\n",
    "            fp <- as.numeric(preplot(fit,p)[2])\n",
    "            dfdp[i,j] <- fp\n",
    "            if (plot==1){\n",
    "                png(filename=sprintf('run%i_p%if%i.png',num,(i-1),(j-1)))\n",
    "                plot(x, y)\n",
    "                lines(fit_orig,col='blue')\n",
    "                abline(v=p,lty=2)\n",
    "                abline(a=-p*fp+f,b=fp,col='red')\n",
    "                dev.off()\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "    dfdp\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dfdp_161 = find_dfdp(data[[1]],seedlist[1],plot=1)\n",
    "dfdp_448 = find_dfdp(data[[2]],seedlist[2],plot=1)\n",
    "dfdp_712 = find_dfdp(data[[3]],seedlist[3],plot=1)\n",
    "dfdp_831 = find_dfdp(data[[4]],seedlist[4],plot=1)\n",
    "dfdp_866 = find_dfdp(data[[5]],seedlist[5],plot=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "setwd('~/Assembly Bias Project/Fisher analysis/Run_062718/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np$savez('dfdp_bolp_19p5_um0_tcub_alpha07.npz',dfdp_161=dfdp_161,dfdp_448=dfdp_448,dfdp_712=dfdp_712,\n",
    "         dfdp_831=dfdp_831,dfdp_866=dfdp_866)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R 3.3",
   "language": "R",
   "name": "ir33"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}