12RootOf2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Welcome to Python in CoCalc!\n",
    "\n",
    "Run a cell via `Shift + Return`. Learn more about [CoCalc Jupyter Notebooks](https://doc.cocalc.com/jupyter.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slice_image(im,M=1,N=1):\n",
    "    tiles = [im[x:x+M,y:y+N].flatten() for x in range(0,im.shape[0],M) for y in range(0,im.shape[1],N)]\n",
    "    return tiles\n",
    "\n",
    "def read_image(fn=\"./images/017.png_wave.jpg\"):\n",
    "    import cv2\n",
    "    return cv2.imread(fn,cv2.COLOR_BGR2RGB) # grayscale\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def constructImageData(Ndim=11,imageFN=\"./data/image.jpg\",nr=10):\n",
    "    import numpy as np\n",
    "    from sklearn.decomposition import PCA\n",
    "    \n",
    "    #import umap\n",
    "    #reducer = umap.UMAP(n_neighbors=15,min_dist=0.1,n_components=2*len(tracks),metric=\"euclidean\",random_state=42)\n",
    "    pca = PCA(n_components=Ndim)\n",
    "    #from sklearn import datasets\n",
    "    import pandas as pd\n",
    "\n",
    "    import numpy as np\n",
    "    from sklearn.impute import SimpleImputer\n",
    "\n",
    "    fn = imageFN\n",
    "    #fn = \"./300/image.jpg\"\n",
    "    #nr = 64\n",
    "    im = read_image(fn)\n",
    "    N = im.shape[1]//nr\n",
    "    mydata = slice_image(im,M=im.shape[0],N=N)\n",
    "    print(len(mydata))\n",
    "    #X = mydata.iloc[:,2:]\n",
    "    Xboston = mydata #mydata.data\n",
    "    X = [x/np.sqrt(np.dot(x,x)) if np.dot(x,x)>0.0 else x for x in Xboston]\n",
    "\n",
    "\n",
    "    from sklearn.decomposition import PCA\n",
    "\n",
    "    nDim =Ndim\n",
    "    pca = PCA(n_components=nDim)\n",
    "    from sklearn.preprocessing import StandardScaler,MinMaxScaler\n",
    "    stdScaler = StandardScaler() \n",
    "    \n",
    "    \n",
    "    X = pca.fit_transform(stdScaler.fit_transform(X))\n",
    "    #y = mydata.iloc[:,1].values\n",
    "    #y = mydata.target\n",
    "    return X\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tsp(data):\n",
    "    import numpy as np\n",
    "    from python_tsp.distances import great_circle_distance_matrix\n",
    "\n",
    "    sources = np.array(data)\n",
    "    distance_matrix = great_circle_distance_matrix(sources)\n",
    "    from python_tsp.distances import tsplib_distance_matrix\n",
    "    from python_tsp.heuristics import solve_tsp_local_search, solve_tsp_simulated_annealing\n",
    "\n",
    "    permutation, distance = solve_tsp_simulated_annealing(distance_matrix)\n",
    "    return permutation, distance\n",
    "\n",
    "\n",
    "def constructData(Ndim=12):\n",
    "    from sklearn.datasets import load_boston\n",
    "    Xboston, y = load_boston(return_X_y=True)\n",
    "    X = [x/np.sqrt(np.dot(x,x)) for x in Xboston]\n",
    "\n",
    "\n",
    "    from sklearn.decomposition import PCA\n",
    "\n",
    "    nDim =Ndim\n",
    "    pca = PCA(n_components=nDim)\n",
    "    from sklearn.preprocessing import StandardScaler,MinMaxScaler\n",
    "    stdScaler = StandardScaler() \n",
    "    \n",
    "    \n",
    "    X = pca.fit_transform(stdScaler.fit_transform(X))\n",
    "    #y = mydata.iloc[:,1].values\n",
    "    #y = mydata.target\n",
    "    return X\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3/512, 3/256, 3/128, 3/64, 3/32, 3/16, 3/8, 3/4, 3/2, 3, 6, 12]\n",
      "12\n",
      "[(0, 0), (1, 1)]\n",
      "(0, 0)\n",
      "(1, 1)\n",
      "written\n"
     ]
    }
   ],
   "source": [
    "#from sage.all import *\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import pandas,sys\n",
    "\n",
    "import statsmodels.api as sm\n",
    "\n",
    "from statsmodels.tsa.api import VAR\n",
    "\n",
    "def parseMidi(fp,part=0):\n",
    "    import os\n",
    "    from music21 import converter\n",
    "    #print(fp)\n",
    "    score = converter.parse(fp,quantizePost=True)\n",
    "    #print(list(score.elements[0].notesAndRests))\n",
    "    score.makeVoices()\n",
    "    #print(dir(score))\n",
    "    #print([e.partAbbreviation for e in score.elements][0])\n",
    "    from music21 import chord\n",
    "    durs = []\n",
    "    ll0 = []\n",
    "    vols = []\n",
    "    isPauses = []\n",
    "    for p in score.elements[part].recurse().notesAndRests: #parts[0].streams()[part].notesAndRests:\n",
    "        #print(p)\n",
    "        if type(p)==chord.Chord:\n",
    "            pitches = sorted([e.pitch.midi-21 for e in p])[0] # todo: think about chords\n",
    "            vol = sorted([e.volume.velocity for e in p])[0]\n",
    "            dur = float(p.duration.quarterLength)\n",
    "            #print(pitches)\n",
    "            \n",
    "            isPause = False\n",
    "        elif (p.name==\"rest\"):\n",
    "            pitches = 89\n",
    "            vol = 64\n",
    "            dur = float(p.duration.quarterLength)\n",
    "            \n",
    "            isPause = True\n",
    "        else:\n",
    "            pitches = p.pitch.midi-21\n",
    "            vol = p.volume.velocity\n",
    "            dur = float(p.duration.quarterLength)\n",
    "            \n",
    "            isPause =  False\n",
    "        if not dur>0 and vol>0:\n",
    "            continue\n",
    "        ll0.append(min(87,pitches))    \n",
    "        durs.append(dur/(12*4.0))\n",
    "        if vol is None or vol == 0:\n",
    "            vol = 64\n",
    "        vols.append(vol*1.0/127.0)\n",
    "        isPauses.append(isPause)\n",
    "        #print(pitches,dur,vol,isPause)\n",
    "            #print(p.name,p.octave,p.duration.quarterLength)\n",
    "    #print(dir(score)) \n",
    "    #print(ll0)\n",
    "    #print(durs)\n",
    "    return ll0,durs,vols,isPauses\n",
    "\n",
    "\n",
    "\n",
    "def dist(k1,k2):\n",
    "    q = getRational(k2-k1)\n",
    "    a,b = q.numerator(),q.denominator()\n",
    "    return sqrt(2*(1-gcd(a,b)^2/(a*b)))\n",
    "\n",
    "def kernPause(a1,a2):\n",
    "    return  1*(a1==a2)\n",
    "\n",
    "def kernPitch(k1,k2):\n",
    "    q = getRational(k2-k1)\n",
    "    a,b = q.numerator(),q.denominator()\n",
    "    return gcd(a,b)**2/(a*b)\n",
    "\n",
    "def kernDuration(k1,k2):\n",
    "    if k1==1.0 or k2==1.0:\n",
    "        return 0.0\n",
    "    return  log(k1)*log(k2)/sqrt(log(k1)**2*log(k2)**2)\n",
    "\n",
    "def kernVolume(v1,v2):\n",
    "    if v1==1.0 or v2==1.0:\n",
    "        return 0.0\n",
    "    return log(v1)*log(v2)/sqrt(log(v1)**2*log(v2)**2)\n",
    "\n",
    "def kern(t1,t2):\n",
    "    pitch1,duration1,volume1,isPause1 = t1\n",
    "    pitch2,duration2,volume2,isPause2 = t2\n",
    "    return 1.0/4*(kernPause(isPause1,isPause2)+kernPitch(pitch1,pitch2)+kernDuration(duration1,duration2)+kernVolume(volume1,volume2))\n",
    "\n",
    "def getRational(k):\n",
    "    alpha = 2**(1/12.0)\n",
    "    x = RDF(alpha**k).n(50)\n",
    "    return x.nearby_rational(max_error=0.01*x)\n",
    "\n",
    "\n",
    "def ngrams(input, n):\n",
    "    output = []\n",
    "    for i in range(len(input)-n+1):\n",
    "        output.append(input[i:i+n])\n",
    "    return output\n",
    "\n",
    "def kernNgram(ngrams1,ngrams2):\n",
    "    return 1.0/len(ngrams1)*sum([ kern(ngrams1[i], ngrams2[i]) for i in range(len(ngrams1))]) \n",
    "\n",
    "\n",
    "#a#lpha = (2**(1/12.0))\n",
    "#qq=[RDF(alpha**k).n(50).nearby_rational(max_error=0.01*RDF(alpha**k).n(50)) for k in range(12)]\n",
    "\n",
    "\n",
    "def getBestLength(fn,Nmaxlen=20,part=1):\n",
    "    pitches, durations, volumes, isPauses = parseMidi(fn,part=part)\n",
    "    ll = []\n",
    "    zz = list(zip(pitches,durations,volumes,isPauses))\n",
    "    #zz2 = list(zip(pitches2,durations2,volumes2,isPauses2))\n",
    "    #print(len((zz)))\n",
    "    #M = matrix([[1]])\n",
    "    Nlen = len(zz)\n",
    "    s = 0\n",
    "    for Nseq in range(1,Nmaxlen):\n",
    "        Z = ngrams([t for t in zz[0:Nmaxlen]],Nseq)\n",
    "        n = len(Z)\n",
    "    #Z.extend(ngrams([t for t in zz2[0:Nmaxlen2]],Nseq))\n",
    "        M = matrix([[kernNgram(t1,t2)-0.5 for t1 in Z] for t2 in Z],ring=RDF)\n",
    "        #print(M)\n",
    "        s=sum(sum(M))/n**2\n",
    "        ll.append((s,Nseq))\n",
    "    return sorted(ll,reverse=True)   \n",
    "      \n",
    "\n",
    "def measureConsonance(pitches,durations,volumes,isPauses,Nmaxlen):\n",
    "    zz = list(zip(pitches,durations,volumes,isPauses))\n",
    "    #zz2 = list(zip(pitches2,durations2,volumes2,isPauses2))\n",
    "    #print(len((zz)))\n",
    "    #M = matrix([[1]])\n",
    "    Nlen = len(zz)\n",
    "    s = 0\n",
    "    for Nseq in range(1,Nmaxlen):\n",
    "        Z = ngrams([t for t in zz[0:Nmaxlen]],Nseq)\n",
    "    #Z.extend(ngrams([t for t in zz2[0:Nmaxlen2]],Nseq))\n",
    "        M = matrix([[kernNgram(t1,t2) for t1 in Z] for t2 in Z],ring=RDF)\n",
    "        #print(M)\n",
    "        s+=sum(sum(M))/(Nmaxlen-Nseq+1)**2.0\n",
    "    return s/(Nmaxlen*1.0)    \n",
    "    \n",
    "def getGramMatrix(pitches,durations,volumes,isPauses,Nmaxlen,Nseq):\n",
    "    zz = list(zip(pitches,durations,volumes,isPauses))\n",
    "    #zz2 = list(zip(pitches2,durations2,volumes2,isPauses2))\n",
    "    print(len((zz)))\n",
    "    #M = matrix([[1]])\n",
    "    Nlen = len(zz)\n",
    "    Z = ngrams([t for t in zz[0:Nmaxlen]],Nseq)\n",
    "    #Z.extend(ngrams([t for t in zz2[0:Nmaxlen2]],Nseq))\n",
    "    M = matrix([[kernNgram(t1,t2) for t1 in Z] for t2 in Z],ring=RDF)\n",
    "    if not M.is_positive_definite():\n",
    "        M+= matrix.identity(len(Z))*0.1\n",
    "    return M\n",
    "    \n",
    "def similarity(fn1,fn2,part=0,Nmaxlen=50):\n",
    "    pitches,durations,volumes,isPauses = parseMidi(fn1,part=part)\n",
    "    pitches2,durations2,volumes2,isPauses2 = parseMidi(fn2,part=part)\n",
    "    \n",
    "    zz = list(zip(pitches,durations,volumes,isPauses))\n",
    "    zz2 = list(zip(pitches2,durations2,volumes2,isPauses2))\n",
    "    #zz2 = list(zip(pitches2,durations2,volumes2,isPauses2))\n",
    "    #print(len((zz)))\n",
    "    #M = matrix([[1]])\n",
    "    Nlen = len(zz)\n",
    "    \n",
    "    \n",
    "    Z = ngrams([t for t in zz[0:Nmaxlen]],Nmaxlen-1)\n",
    "    Z2 = ngrams([t for t in zz2[0:Nmaxlen]],Nmaxlen-1)\n",
    "    \n",
    "    return kernNgram(Z[0],Z2[0])\n",
    "    \n",
    "def getVarModel(fn,voice=1,Nmaxlen=100,Ndim=3,Nseq=8,maxlags=6,nCopy=4):\n",
    "    pitches,durations,volumes,isPauses = parseMidi(fn,part=voice)\n",
    "    \n",
    "    #print(\"consonance: \",measureConsonance(pitches,durations,volumes,isPauses,40))\n",
    "    #pitches2,durations2,volumes2,isPauses2 = parseMidi(fn2,part=voice)\n",
    "    \n",
    "    zz = list(zip(pitches,durations,volumes,isPauses))\n",
    "    #zz2 = list(zip(pitches2,durations2,volumes2,isPauses2))\n",
    "    #print(len((zz)))\n",
    "    #M = matrix([[1]])\n",
    "    Nlen = len(zz)\n",
    "    \n",
    "    ZZ = []\n",
    "\n",
    "    ZZ.extend(zz[0:Nmaxlen])\n",
    "    #ZZ.extend(zz[(Nlen//2-Nmaxlen):(Nlen//2+Nmaxlen)])\n",
    "    #ZZ.extend(zz[-Nmaxlen:])\n",
    "    \n",
    "    Z = ngrams([t for t in ZZ],Nseq)\n",
    "    #Z.extend(ngrams([t for t in zz2[0:Nmaxlen2]],Nseq))\n",
    "    M = matrix([[kernNgram(t1,t2) for t1 in Z] for t2 in Z],ring=RDF)\n",
    "    if not M.is_positive_definite():\n",
    "        M+= matrix.identity(len(Z))*0.1\n",
    "\n",
    "#print(qq)\n",
    "#M = matrix([[kern(k1,k2) for k1 in range(88)] for k2 in range(88)])\n",
    "    #print(M.str())\n",
    "#print(M.is_positive_definite())\n",
    "\n",
    "    #print(M)\n",
    "    CC = M.cholesky()\n",
    "#Ch = M.cholesky().rows()\n",
    "\n",
    "    from sklearn.decomposition import PCA\n",
    "\n",
    "    nDim =Ndim\n",
    "    pca = PCA(n_components=nDim)\n",
    "    from sklearn.preprocessing import StandardScaler,MinMaxScaler\n",
    "    stdScaler = StandardScaler()\n",
    "    Ch = pca.fit_transform(stdScaler.fit_transform(CC))\n",
    "\n",
    "    model = VAR(Ch)\n",
    "\n",
    "    results = model.fit(maxlags=maxlags,trend=\"nc\",ic=\"aic\")\n",
    "    #results = model.fit(maxlags,trend=\"nc\")\n",
    "    #print(results.summary())\n",
    "\n",
    "    lag_order = results.k_ar\n",
    "\n",
    "    preds = []\n",
    "    X = [x for x in Ch]\n",
    "    return results, X,Ch,Z,CC\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def getCholesky(fn,voice=1,Nmaxlen=100,Ndim=3,Nseq=8,maxlags=6,nCopy=4):\n",
    "    pitches,durations,volumes,isPauses = parseMidi(fn,part=voice)\n",
    "    \n",
    "    #print(\"consonance: \",measureConsonance(pitches,durations,volumes,isPauses,40))\n",
    "    #pitches2,durations2,volumes2,isPauses2 = parseMidi(fn2,part=voice)\n",
    "    \n",
    "    zz = list(zip(pitches,durations,volumes,isPauses))\n",
    "    #zz2 = list(zip(pitches2,durations2,volumes2,isPauses2))\n",
    "    print(\"len(zz) = \",len((zz)))\n",
    "    #M = matrix([[1]])\n",
    "    Nlen = len(zz)\n",
    "    \n",
    "    ZZ = []\n",
    "\n",
    "    ZZ.extend(zz[0:Nmaxlen])\n",
    "    #ZZ.extend(zz[(Nlen//2-Nmaxlen):(Nlen//2+Nmaxlen)])\n",
    "    #ZZ.extend(zz[-Nmaxlen:])\n",
    "    \n",
    "    Z = ngrams([t for t in ZZ],Nseq)\n",
    "    #Z.extend(ngrams([t for t in zz2[0:Nmaxlen2]],Nseq))\n",
    "    M = matrix([[kernNgram(t1,t2) for t1 in Z] for t2 in Z],ring=RDF)\n",
    "    if not M.is_positive_definite():\n",
    "        M+= matrix.identity(len(Z))*0.01\n",
    "\n",
    "#print(qq)\n",
    "#M = matrix([[kern(k1,k2) for k1 in range(88)] for k2 in range(88)])\n",
    "    #print(M.str())\n",
    "#print(M.is_positive_definite())\n",
    "\n",
    "    #print(M)\n",
    "    CC = M.cholesky()\n",
    "#Ch = M.cholesky().rows()\n",
    "\n",
    "    from sklearn.decomposition import PCA\n",
    "\n",
    "    nDim =Ndim\n",
    "    pca = PCA(n_components=nDim)\n",
    "    from sklearn.preprocessing import StandardScaler,MinMaxScaler\n",
    "    stdScaler = StandardScaler()\n",
    "    Ch = pca.fit_transform(stdScaler.fit_transform(CC))\n",
    "\n",
    "    #model = VAR(Ch)\n",
    "\n",
    "    #results = model.fit(maxlags=maxlags,trend=\"nc\",ic=\"aic\")\n",
    "    #results = model.fit(maxlags,trend=\"nc\")\n",
    "    #print(results.summary())\n",
    "\n",
    "    #lag_order = results.k_ar\n",
    "\n",
    "    preds = []\n",
    "    X = [x for x in Ch]\n",
    "    return  X,Ch,Z,CC\n",
    "\n",
    "\n",
    "def kk(a,b):\n",
    "    return gcd(a,b)^2/(a*b)\n",
    "\n",
    "def getData(Ndim,numbers=range(1,100),kern=kk):\n",
    "    M = matrix([[kern(a,b) for a in numbers] for b in numbers],ring=RDF)\n",
    "    if not M.is_positive_definite():\n",
    "        M += 0.1*matrix.identity(len(numbers))\n",
    "    CC = M.cholesky()\n",
    "    \n",
    "    from sklearn.decomposition import PCA\n",
    "    nDim =Ndim\n",
    "    pca = PCA(n_components=nDim)\n",
    "    from sklearn.preprocessing import StandardScaler,MinMaxScaler\n",
    "    stdScaler = StandardScaler()\n",
    "    Ch = pca.fit_transform(stdScaler.fit_transform(CC))\n",
    "    X = [x for x in Ch]\n",
    "    return X\n",
    "    \n",
    "def forecastData(fn,voice=1,Nmaxlen=100,Ndim=3,Nseq=10,maxlags=6,nCopy=4,numbers=range(1,100),kern=kk,useThisDataX=None):\n",
    "    X1,Ch1,Z1,CC1 = getCholesky(fn,voice[0],Nmaxlen,Ndim,Nseq,maxlags,nCopy)\n",
    "    #result2,X2,Ch2,Z2,CC2 = getVarModel(fn2,voice[1],Nmaxlen,Ndim,Nseq,maxlags,nCopy)\n",
    "    if useThisDataX is None:\n",
    "        X1 = getData(Ndim,numbers,kern)\n",
    "        Nforecast = len(numbers)\n",
    "    else:\n",
    "        X1 = useThisDataX\n",
    "        Nforecast = len(X1)\n",
    "    #lag_order = result1.k_ar\n",
    "    preds = []\n",
    "    \n",
    "    \n",
    "    \n",
    "    for k in range(Nforecast):\n",
    "        preds.append(X1[k])\n",
    "\n",
    "    from sklearn.neighbors import NearestNeighbors\n",
    "\n",
    "    nbrs = NearestNeighbors(n_neighbors=1, algorithm='ball_tree').fit(Ch1)\n",
    "\n",
    "    def findBestMatch(nbrs,new_row):\n",
    "        distances,indices = nbrs.kneighbors([np.array(new_row)])\n",
    "        return indices[0][0]\n",
    "\n",
    "\n",
    "    inds = []\n",
    "    for x in preds:\n",
    "        i = findBestMatch(nbrs,x)\n",
    "        \n",
    "        lz = len(Z1[i])//2\n",
    "        \n",
    "        inds.extend(Z1[i][lz-nCopy:lz+nCopy+1])\n",
    "\n",
    "    return(inds)\n",
    "\n",
    "    \n",
    "    \n",
    "def forecastPart(fn,fn2,voice=1,Nforecast=20,Nmaxlen=100,Ndim=3,Nseq=8,maxlags=6,nCopy=4):\n",
    "    result1,X1,Ch1,Z1,CC1 = getVarModel(fn,voice[0],Nmaxlen,Ndim,Nseq,maxlags,nCopy)\n",
    "    result2,X2,Ch2,Z2,CC2 = getVarModel(fn2,voice[1],Nmaxlen,Ndim,Nseq,maxlags,nCopy)\n",
    "    #X1 = getData(Nforecast,Ndim)\n",
    "    lag_order = result2.k_ar\n",
    "    preds = []\n",
    "    \n",
    "    \n",
    "    for k in range(Nforecast):\n",
    "        a2 = 0.9\n",
    "        pred = result2.forecast(np.array((X1[k:(k+lag_order)])), 1)\n",
    "        preds.append(pred[0])\n",
    "        X1.append(pred[0])\n",
    "\n",
    "    from sklearn.neighbors import NearestNeighbors\n",
    "\n",
    "    nbrs = NearestNeighbors(n_neighbors=1, algorithm='ball_tree').fit(Ch1)\n",
    "\n",
    "    def findBestMatch(nbrs,new_row):\n",
    "        distances,indices = nbrs.kneighbors([np.array(new_row)])\n",
    "        return indices[0][0]\n",
    "\n",
    "\n",
    "    inds = []\n",
    "    for x in preds:\n",
    "        i = findBestMatch(nbrs,x)       \n",
    "        lz = len(Z1[i])//2\n",
    "        inds.extend(Z1[i][lz-nCopy:lz+nCopy+1])\n",
    "\n",
    "    return(inds)\n",
    "\n",
    "durationslist = [[sum([(2**(n-i)) for i in range(d+1)]) for n in range(-8,3+1)] for d in range(1,2)]\n",
    "notevalues = []\n",
    "for i in range(len(durationslist)):\n",
    "    notevalues.extend(durationslist[i])\n",
    "    \n",
    "notevalues = sorted(notevalues)    \n",
    "print(notevalues)\n",
    "print(len(notevalues))\n",
    "\n",
    "def findNearestDuration(duration):\n",
    "    return sorted([(abs(duration-nv),nv) for nv in notevalues])[0][1]\n",
    "\n",
    "def writePitches(fn,inds,tempo=82,instrument=[0,0]):\n",
    "    from MidiFile import MIDIFile\n",
    "\n",
    "    track    = 0\n",
    "    channel  = 0\n",
    "    time     = 0   # In beats\n",
    "    duration = 1   # In beats\n",
    "    tempo    = tempo # In BPM\n",
    "    volume   = 116 # 0-127, as per the MIDI standard\n",
    "\n",
    "    ni = len(inds)\n",
    "    MyMIDI = MIDIFile(ni,adjust_origin=False) # One track, defaults to format 1 (tempo track\n",
    "                     # automatically created)\n",
    "    MyMIDI.addTempo(track,time, tempo)\n",
    "\n",
    "\n",
    "    for k in range(ni):\n",
    "        MyMIDI.addProgramChange(k,k,0,instrument[k])\n",
    "\n",
    "\n",
    "    times = ni*[0]\n",
    "    for k in range(len(inds)):\n",
    "        channel = k\n",
    "        track = k\n",
    "        for i in range(len(inds[k])):\n",
    "            pitch,duration,volume,isPause = inds[k][i]\n",
    "            track = k\n",
    "            channel = k\n",
    "            duration = duration*12*4 #findNearestDuration(duration*12*4)\n",
    "            #print(k,pitch,times[k],duration,100)\n",
    "            if not isPause: #rest\n",
    "                #print(volumes[i])\n",
    "                # because of median:\n",
    "                pitch = int(floor(pitch))\n",
    "                MyMIDI.addNote(track, channel, pitch+21, float(times[k]) , float(duration), int(127*volume))#*(ni-k+10.0)/(ni+10.0))\n",
    "            times[k] += duration*1.0    \n",
    "       \n",
    "    with open(fn, \"wb\") as output_file:\n",
    "        MyMIDI.writeFile(output_file)\n",
    "    print(\"written\")    \n",
    " \n",
    "\n",
    "\n",
    "\n",
    "def getMostSimilar(fn,fns,parts = [0,1]):\n",
    "    return sorted([(sum([similarity(fn,f,part=part) for part in parts]),f) for f in fns if f!=fn],reverse=True)[0][1]\n",
    "\n",
    "\n",
    "\n",
    "def midiWithData(fn,pps,instruments=[0]*3,numbers=range(1,100),kern=kk,useThisDataX=None,Ncopys=[0]):\n",
    "    \n",
    "    iinds = []\n",
    "    for i in range(len(pps)):\n",
    "        p = pps[i]\n",
    "        print(p)\n",
    "        inds = forecastData(fn,voice=p,Nmaxlen=200,\n",
    "                            Ndim=11,Nseq=19,maxlags=6,nCopy=Ncopys[i],\n",
    "                            numbers=numbers[i],kern=kern,useThisDataX=useThisDataX)\n",
    "        iinds.append(inds)\n",
    "    fn = fn+\".mix.mid\"    \n",
    "    writePitches(fn,iinds,tempo=80,instrument=instruments)\n",
    "\n",
    "def rewrite(str,rules):\n",
    "    s = \"\"\n",
    "    for x in str:\n",
    "        s += rules[x]\n",
    "    return s\n",
    "\n",
    "def iterRewrite(str,rules,n=1):\n",
    "    s = str\n",
    "    for x in range(n):\n",
    "        s = rewrite(s,rules)\n",
    "    return s    \n",
    "\n",
    "    \n",
    "def mixMidis(fn1,fn2,instruments=[0]*9,parts=list(zip(list(range(2)),[0,1])), nCopy=[4,0]):\n",
    "    #pps = [0,1,2,3,4,5]\n",
    "    pps = parts \n",
    "    print(pps)\n",
    "    iinds = []\n",
    "    for i in range(len(pps)):\n",
    "        p = pps[i]\n",
    "        print(p)\n",
    "        inds = forecastPart(fn1,fn2,voice=p,Nforecast=200,Nmaxlen=200,Ndim=7,Nseq=7,maxlags=4,nCopy=nCopy[i])\n",
    "    #print(measureConsonance(pitches,durations,volumes,isPauses,Nmaxlen=23))\n",
    "    #print(inds)\n",
    "        iinds.append(inds)\n",
    "    fn1 = fn1.split(\"/\")[-1]\n",
    "    fn2 = fn2.split(\"/\")[-1]\n",
    "    fn = \"./midi/\"+fn1+\"_\"+fn2+\".mix.mid\"    \n",
    "    writePitches(fn,iinds,tempo=80,instrument=instruments) #[0]*len(iinds))    \n",
    "\n",
    "#from constants import INSTRUMENT_MAP as im, DRUM_MAP as dm\n",
    "\n",
    "instruments = []\n",
    "\n",
    "def kk(a,b):\n",
    "    return gcd(a,b)**2/(a*b)\n",
    "\n",
    "def kk2(a,b):\n",
    "    return 2*gcd(a,b)/(a+b)\n",
    "\n",
    "def kk3(a,b):\n",
    "    return 2*gcd(a,b)**3/(a*b*(a+b))\n",
    "\n",
    "\n",
    "def mainNumbersData():\n",
    "    pass\n",
    "\n",
    "def mainImageData(fn):\n",
    "    imageFN = \"./data/two.jpg\"\n",
    "    div = divisors(read_image(imageFN).shape[1]) # divisors of the length in y-direction of the image\n",
    "    nr = sorted([ (abs(d-100),d) for d in div])[0][1] # the closest divisor to the number 50\n",
    "    print(nr)\n",
    "    data = constructImageData(Ndim=19,imageFN=imageFN,nr=nr)\n",
    "    useThisDataX = data\n",
    "    pps = list(zip([0,1],[0,0]))\n",
    "    kk3 = None\n",
    "    midiWithData(fn,pps,instruments=2*[0],numbers=numbers,kern=kk3,useThisDataX=useThisDataX,Ncopys=[3,1])\n",
    "\n",
    "\n",
    "def mainTSPData(fn):\n",
    "    data = constructData(Ndim=11) \n",
    "    permutation,distance = tsp(data)\n",
    "    data = [data[permutation[k]] for k in range(len(data))]\n",
    "    useThisDataX = data\n",
    "    pps = list(zip([0,1],[0,0]))\n",
    "    kk3 = None\n",
    "    midiWithData(fn,pps,instruments=2*[0],numbers=numbers,kern=kk3,useThisDataX=useThisDataX,Ncopys=[3,1])\n",
    "\n",
    "\n",
    "def mainMixMidis(fn1,fn2):\n",
    "    mixMidis(fn1,fn2,instruments=2*[0],parts=list(zip([0,1],[0,1])),nCopy=[3,2])\n",
    "\n",
    "\n",
    "def main():\n",
    "    import os\n",
    "    fn1 = sorted([os.path.join(\"./midi/\",fn) for fn in os.listdir(\"./midi/\") if fn.endswith(\".mid\") \n",
    "              and \"clair_de_lune_theme\" in fn.lower() and not \"mix\" in fn.lower()])[0]\n",
    "    \n",
    "    fn2 = sorted([os.path.join(\"./midi/\",fn) for fn in os.listdir(\"./midi/\") if fn.endswith(\".mid\") \n",
    "              and \"swan_theme\" in fn.lower() and not \"mix\" in fn.lower() and not \"mix\" in fn.lower()])[0]\n",
    "    \n",
    "    #mainTSPData(fn1)\n",
    "    #mainImageData(fn1)\n",
    "    mainMixMidis(fn1,fn2)\n",
    "    \n",
    "\n",
    "main()    \n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]]\n"
     ]
    }
   ],
   "source": [
    "print(list(Permutations([0,1,2])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 14, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 14, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 14, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 14, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 14, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 13, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 14, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 13, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 14, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 13, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 13, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 14, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 14, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 14, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26, 25, 24, 22, 19, 16, 13, 11, 8, 5, 3, 2, 1, 1, 1, 2, 3, 5, 8, 11, 13, 16, 19, 22, 24, 25, 26, 27, 26]\n"
     ]
    }
   ],
   "source": [
    "print([int(floor(13*sin(x*pi.n()/(2*7)).n()))+14 for x in range(1,401)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 0,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.0",
   "language": "sage",
   "name": "sagemath"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}