mixed-use decomposition.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import csv\n",
    "from scipy import optimize as op\n",
    "import math\n",
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#input basedcurves and observed TASs\n",
<<<<<<< HEAD
    "path=\"E:/\"\n",
=======
    "path=\"E:/research/\"\n",
>>>>>>> a2a76978808ca32aa860d3c9c879bf81d420025c
    "\n",
    "observedpath=path+\"observed.csv\"\n",
    "basecurvepath=path+\"basecurve.csv\"\n",
    "\n",
    "observed=np.array(list(csv.reader(open(observedpath,'r'))))\n",
    "basecurve=np.array(list(csv.reader(open(basecurvepath,'r'))))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cosine similarity\n",
    "def calcos(A,B):\n",
    "    N=len(A)\n",
    "    suma2=0\n",
    "    sumb2=0\n",
    "    sumab=0\n",
    "    for i in range(N):\n",
    "        suma2+=float(A[i])*float(A[i])\n",
    "        sumb2+=float(B[i])*float(B[i])\n",
    "        sumab+=float(A[i])*float(B[i])\n",
    "    if(suma2==0 or sumb2==0):\n",
    "        return -2\n",
    "    r=sumab/(math.sqrt(suma2)*math.sqrt(sumb2))\n",
    "    return r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#optimize\n",
    "#bases and tmpcurve are global variables\n",
    "def similarity(X):\n",
    "    L=len(tmpcurve)\n",
    "    N=len(X)\n",
    "    newcurve=np.zeros(L)\n",
    "    for i in range(N):\n",
    "        for j in range(L):\n",
    "            newcurve[j]+=float(X[i])*float(bases[i][j])\n",
    "    result=calcos(newcurve,tmpcurve)\n",
    "    return 1-result\n",
    "\n",
    "def constraint(X):\n",
    "    return sum(X)-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#other results\n",
    "#mixing index\n",
    "def calmix(weights):\n",
    "    result=0\n",
    "    for i in range(len(weights)):\n",
    "        if(float(weights[i])==0):\n",
    "            continue\n",
    "        result+=(-1)*float(weights[i])*math.log(float(weights[i]))\n",
    "    return result\n",
    "\n",
    "#dynamic activity proportion\n",
    "def caldynamic(weights,bases):\n",
    "    L=len(bases[0])\n",
    "    M=len(weights)\n",
    "    total=np.zeros(L)\n",
    "    for i in range(M):\n",
    "        for j in range(L):\n",
    "            total[j]+=float(weights[i])*float(bases[i][j])\n",
    "    result=[]\n",
    "    for j in range(L):\n",
    "        tmp=[]\n",
    "        for i in range(M):\n",
    "            tmpro=float(weights[i])*float(bases[i][j])/total[j]\n",
    "            tmp.append(tmpro)\n",
    "        result.append(tmp)\n",
    "    return np.array(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  0.00000000e+00   9.29515522e-01   0.00000000e+00   7.04844780e-02\n",
      "   6.26402248e-17] 0.254890342183164\n",
      "-------------------------------\n",
      "[  0.00000000e+00   1.00000000e+00   4.67659656e-17   6.93889390e-18\n",
      "   7.88823115e-17] 4.957459789324139e-15\n",
      "-------------------------------\n",
      "[  3.70999641e-01   6.29000359e-01   0.00000000e+00   8.67361738e-17\n",
      "   1.24303779e-16] 0.6594855639200745\n",
      "-------------------------------\n",
      "[  0.00000000e+00   1.00000000e+00   0.00000000e+00   0.00000000e+00\n",
      "   4.86823289e-17] 1.8285674506587357e-15\n",
      "-------------------------------\n",
      "[  0.00000000e+00   6.26826109e-01   8.85630418e-18   1.06278565e-01\n",
      "   2.66895325e-01] 0.8835672430439206\n",
      "-------------------------------\n",
      "[ 0.          0.41317815  0.06223192  0.          0.52458993] 0.8764425530629117\n",
      "-------------------------------\n",
      "[  0.00000000e+00   6.84084203e-01   3.15915797e-01   0.00000000e+00\n",
      "   1.60426812e-18] 0.6237524842820085\n",
      "-------------------------------\n",
      "[  4.16333634e-17   0.00000000e+00   6.93889390e-17   1.00000000e+00\n",
      "   0.00000000e+00] 4.152052458276665e-15\n",
      "-------------------------------\n",
      "[  8.22339438e-17   0.00000000e+00   0.00000000e+00   8.91242774e-01\n",
      "   1.08757226e-01] 0.34390910474889425\n",
      "-------------------------------\n",
      "[  0.00000000e+00   4.10866451e-17   5.52324478e-01   0.00000000e+00\n",
      "   4.47675522e-01] 0.6876614400446317\n",
      "-------------------------------\n",
      "[ 0.04156204  0.47522115  0.21024585  0.27297096  0.        ] 1.1680399964400574\n",
      "[[ 0.04245671  0.46447167  0.18526546  0.30780616  0.        ]\n",
      " [ 0.04264815  0.44606984  0.16847845  0.34280357  0.        ]\n",
      " [ 0.04408979  0.43450062  0.17618003  0.34522956  0.        ]\n",
      " [ 0.04422347  0.42959656  0.17620054  0.34997943  0.        ]\n",
      " [ 0.04596486  0.49694974  0.13379766  0.32328774  0.        ]\n",
      " [ 0.03915659  0.64993721  0.0855264   0.2253798   0.        ]\n",
      " [ 0.03042947  0.6732993   0.07360082  0.22267041  0.        ]\n",
      " [ 0.02735949  0.64795297  0.07348081  0.25120673  0.        ]\n",
      " [ 0.03114217  0.58121412  0.07369642  0.31394729  0.        ]\n",
      " [ 0.04270063  0.50830058  0.07940598  0.36959281  0.        ]\n",
      " [ 0.05260774  0.50329418  0.09392374  0.35017433  0.        ]\n",
      " [ 0.05208477  0.45065649  0.19681683  0.30044191  0.        ]\n",
      " [ 0.04186192  0.37452134  0.33565823  0.24795851  0.        ]\n",
      " [ 0.04730588  0.4158789   0.2570783   0.27973692  0.        ]\n",
      " [ 0.05450633  0.46125983  0.1867764   0.29745744  0.        ]\n",
      " [ 0.05621754  0.47532611  0.16808514  0.30037121  0.        ]\n",
      " [ 0.0521929   0.49660621  0.17537024  0.27583065  0.        ]\n",
      " [ 0.03858005  0.48755419  0.23336471  0.24050106  0.        ]\n",
      " [ 0.02921194  0.43574494  0.32511397  0.20992916  0.        ]\n",
      " [ 0.03028213  0.39881967  0.37476831  0.19612989  0.        ]\n",
      " [ 0.0355398   0.4189164   0.33497067  0.21057312  0.        ]\n",
      " [ 0.0417391   0.46624316  0.24938602  0.24263173  0.        ]\n",
      " [ 0.04217973  0.48950792  0.1960059   0.27230645  0.        ]\n",
      " [ 0.04278663  0.46505836  0.17526576  0.31688924  0.        ]]\n",
      "-------------------------------\n",
      "[ 0.          0.58690758  0.18520785  0.22788457  0.        ] 0.9620897647005595\n",
      "-------------------------------\n",
      "[  7.75791923e-18   3.91715058e-01   2.12269649e-01   2.89051932e-01\n",
      "   1.06963362e-01] 1.2939680681603245\n",
      "-------------------------------\n"
     ]
    }
   ],
   "source": [
    "#decomposition\n",
    "decomposed=np.array(observed)\n",
    "bases=np.array(basecurve)\n",
    "N=len(decomposed)\n",
    "M=len(bases)\n",
    "\n",
    "#>=0,<=1\n",
    "b=(0.0,1.0)\n",
    "bnds=[]\n",
    "#initialization weights\n",
    "x0=[]\n",
    "for i in range(M):\n",
    "    tmp=(0.0,1.0)\n",
    "    bnds.append(tmp)\n",
    "    x0.append(1.0/M)\n",
    "bnds=tuple(bnds)\n",
    "\n",
    "\n",
    "cons=[{'type' : 'eq' , 'fun' : constraint}]\n",
    "\n",
    "for j in range(N):\n",
    "    tmpcurve=copy.deepcopy(decomposed[j])\n",
    "    sol=op.minimize(similarity,x0,method='SLSQP',bounds=bnds,constraints=cons)\n",
    "    composition=np.array(sol['x'])\n",
    "    dynamicpro=caldynamic(composition,bases)\n",
    "    mixingindex=calmix(composition)\n",
    "    print(composition,mixingindex)\n",
    "    if(j==10):\n",
    "        print(dynamicpro)\n",
    "    print(\"-------------------------------\")\n",
    "\n"
   ]
  },
  {
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}