bezier_curve.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 310,
   "id": "1cf0f5fd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1/256, 1/128, 1/64, 1/32, 1/16, 1/8, 1/4, 1/2, 1, 3/512, 3/256, 3/128, 3/64, 3/32, 3/16, 3/8, 3/4, 3/2]\n"
     ]
    }
   ],
   "source": [
    "from sage.all import *\n",
    "import numpy as np\n",
    "import random\n",
    "\n",
    "def kernPause(a1,a2):\n",
    "    return  1*(a1==a2)\n",
    "\n",
    "def kernDuration(d1,d2):\n",
    "    return min(d1,d2)/max(d1,d2)\n",
    "\n",
    "def kernVolume(v1,v2):\n",
    "    return min(v1,v2)/max(v1,v2)\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 kernNote(n1,n2):\n",
    "    p1,d1,v1,r1 = n1\n",
    "    p2,d2,v2,r2 = n2\n",
    "    #print(n1,n2)\n",
    "    return 1/4*(kernPitch(p1,p2)+min(d1,d2)/max(d1,d2)+min(v1,v2)/max(v1,v2)+kernPause(r1,r2))\n",
    "\n",
    "import portion as PP\n",
    "\n",
    "def muInterval(i):\n",
    "    if i.empty:\n",
    "        return 0\n",
    "    return i.upper-i.lower\n",
    "\n",
    "def jaccard(i1,i2):\n",
    "    return muInterval(i1.intersection(i2))/muInterval(i1.union(i2))\n",
    "\n",
    "def kernJacc(interval1,interval2):\n",
    "    min1,max1 = interval1\n",
    "    min2,max2 = interval2\n",
    "    eps = 1.0/128.0\n",
    "    X = PP.closed(min1-eps,max1+eps)\n",
    "    Y = PP.closed(min2-eps,max2+eps)\n",
    "    return jaccard(X,Y)\n",
    "\n",
    "def kernInt(i1,i2):\n",
    "    min1, max1,note1 = i1 \n",
    "    min2, max2,note2 = i2\n",
    "    return kernNote(note1,note2)*kernJacc((min1,max1),(min2,max2))\n",
    "\n",
    "\n",
    "def distNote(kernNote):\n",
    "    import numpy as np\n",
    "    return (lambda n1,n2: np.sqrt(2.0-2*kernNote(n1,n2)))\n",
    "\n",
    "def kernVolume(v1,v2):\n",
    "    #return kernJacc(v1,v2)\n",
    "    return min(v1,v2) #/max(v1,v2)\n",
    "\n",
    "def kernDuration(d1,d2):\n",
    "    return min(d1,d2)\n",
    "\n",
    "def create_random_graph(n):\n",
    "    G = Graph(loops=False)\n",
    "    G.add_vertex(1-1)\n",
    "    for k in range(2,n+1):\n",
    "        vert = [v for v in G.vertices()]\n",
    "        G.add_vertex(k-1)\n",
    "        for v in vert:\n",
    "            prob = 1.0/v*k/sigma(k)\n",
    "            p = randint(1,100)/100.0\n",
    "            #print(prob,p)\n",
    "            if p <= prob and k%(v+1)==0:\n",
    "                G.add_edge(v,k-1)\n",
    "    return G\n",
    "\n",
    "def create_divisor_graph(n):\n",
    "    G = DiGraph(loops=True)\n",
    "    divs = divisors(n)\n",
    "    for d in divs:\n",
    "        G.add_vertex(divs.index(d))\n",
    "    for u in divs:\n",
    "        for v in divs:\n",
    "            if is_prime(u%v):\n",
    "                G.add_edge(divs.index(u),divs.index(v))\n",
    "    return G\n",
    "\n",
    "def kernAdd(t1,t2,alphaPitch=0.25):\n",
    "    pitch1,volume1,isPause1 = t1\n",
    "    pitch2,volume2,isPause2 = t2\n",
    "    #return 1.0/3*(1-alphaPitch)*kernPause(isPause1,isPause2)+alphaPitch*kernPitch(pitch1,pitch2)+1.0/3*(1-alphaPitch)*kernDuration(duration1,duration2)+1.0/3*(1-alphaPitch)*kernVolume(volume1,volume2)\n",
    "    apa = alphaPitch[\"pause\"]\n",
    "    api = alphaPitch[\"pitch\"]\n",
    "    avo = alphaPitch[\"volume\"]\n",
    "    #return kernPause(isPause1,isPause2)*kernPitch(pitch1,pitch2)*kernVolume(volume1,volume2)\n",
    "\n",
    "    if np.abs(apa+api+avo-1)<10**-5:\n",
    "        return apa*kernPause(isPause1,isPause2)+api*kernPitch(pitch1,pitch2)+avo*kernVolume(volume1,volume2)\n",
    "    else:\n",
    "        return None\n",
    "\n",
    "def kern0(zz0,alphaPitch={\"pitch\":1,\"volume\":2,\"pause\":3}):\n",
    "    return lambda t1,t2: kernAdd(zz0[int(t1[0])],zz0[int(t2[0])],alphaPitch)\n",
    "\n",
    "def kern(alphaPitch={\"pitch\":1,\"volume\":2,\"pause\":3}):\n",
    "    return lambda t1,t2: kernAdd(t1,t2,alphaPitch)\n",
    "\n",
    "\n",
    "\n",
    "def distKern1(x,y,alphaPitch={\"pitch\":1,\"volume\":2,\"pause\":3}):\n",
    "    #print(alphaPitch)\n",
    "    return np.sqrt(2-2*kern(alphaPitch)(x,y))\n",
    "\n",
    "def distKern(kern):\n",
    "    return lambda a,b : np.sqrt(kern(a,a)+kern(b,b)-2*kern(a,b))\n",
    "\n",
    "import music21 as m21\n",
    "from itertools import product\n",
    "\n",
    "durlist = [[sum([(power(2,(n-i))) for i in range(d+1)]) for n in range(-8,3+1)] for d in range(2)]\n",
    "durationslist = []\n",
    "for dl in durlist:\n",
    "    durationslist.extend([x for x in dl if x<=1.5])\n",
    "print(durationslist)\n",
    "\n",
    "def findNearestDuration(duration,durationslist):\n",
    "    return sorted([(abs(duration-nv),nv) for nv in durationslist])[0][1]\n",
    "\n",
    "def parse_file(xml):\n",
    "    xml_data = m21.converter.parse(xml)\n",
    "    score = []\n",
    "    for part in xml_data.parts:\n",
    "        parts = []\n",
    "        for note in part.recurse().notesAndRests:\n",
    "            if note.isRest:\n",
    "                start = note.offset\n",
    "                duration = float(note.quarterLength)/4.0\n",
    "                vol = 32 #note.volume.velocity\n",
    "                pitch = 60\n",
    "                parts.append([pitch,findNearestDuration(duration,durationslist),vol,True])\n",
    "            elif note.isChord:\n",
    "                note = [n for n in note][1]\n",
    "                start = note.offset\n",
    "                duration = float(note.quarterLength)/4.0\n",
    "                pitch = note.pitch.midi\n",
    "                #print(pitch,duration,note.volume)\n",
    "                vol = note.volume.velocity\n",
    "                if vol is None:\n",
    "                    vol = int(note.volume.realized * 127)\n",
    "                parts.append([pitch,findNearestDuration(duration,durationslist),vol,False])    \n",
    "            else:\n",
    "                #print(note)\n",
    "                start = note.offset\n",
    "                duration = float(note.quarterLength)/4.0\n",
    "                pitch = note.pitch.midi\n",
    "                #print(pitch,duration,note.volume)\n",
    "                vol = note.volume.velocity\n",
    "                if vol is None:\n",
    "                    vol = int(note.volume.realized * 127)\n",
    "                parts.append([pitch,findNearestDuration(duration,durationslist),vol,False])    \n",
    "        score.append(parts)        \n",
    "    return score\n",
    "\n",
    "\n",
    "def writePitches(fn,inds,tempo=82,instrument=[0,0],add21=True,start_at= [0,0],durationsInQuarterNotes=False):\n",
    "    from MidiFile import MIDIFile\n",
    "\n",
    "    track    = 0\n",
    "    channel  = 0\n",
    "    time     = 0   # In beats\n",
    "    duration = 1   # In beats # 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 = start_at\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",
    "            #print(pitch,duration,volume,isPause)\n",
    "            track = k\n",
    "            channel = k\n",
    "            if not durationsInQuarterNotes:\n",
    "                duration = 4*duration#*maxDurations[k] #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",
    "                if add21:\n",
    "                    pitch += 21\n",
    "                #print(pitch,times[k],duration,volume,isPause)    \n",
    "                MyMIDI.addNote(track, channel, int(pitch), float(times[k]) , float(duration), int(volume))\n",
    "                times[k] += duration*1.0  \n",
    "            else:\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",
    "def writeIntervals(fn,inds,tempo=82,instrument=[0,0],add21=True,start_at= [0,0],durationsInQuarterNotes=False):\n",
    "    from MidiFile import MIDIFile\n",
    "\n",
    "    track    = 0\n",
    "    channel  = 0\n",
    "    time     = 0   # In beats\n",
    "    duration = 1   # In beats # 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 = start_at\n",
    "    for k in range(len(inds)):\n",
    "        channel = k\n",
    "        track = k\n",
    "        for i in range(len(inds[k])):\n",
    "            start,end,note = inds[k][i]\n",
    "            pitch,duration,volume,isPause = note\n",
    "            #print(pitch,duration,volume,isPause)\n",
    "            track = k\n",
    "            channel = k\n",
    "            if not durationsInQuarterNotes:\n",
    "                duration = 4*duration#*maxDurations[k] #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",
    "                if add21:\n",
    "                    pitch += 21\n",
    "                #print(pitch,times[k],duration,volume,isPause)    \n",
    "                MyMIDI.addNote(track, channel, int(pitch), float(start) , float(duration), int(volume))\n",
    "                times[k] += duration*1.0  \n",
    "            else:\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",
    "def run_length_row(row_of_01):\n",
    "    from itertools import groupby\n",
    "    ar = row_of_01\n",
    "    return [(k, sum(1 for i in g)) for k,g in groupby(ar)]\n",
    "\n",
    "def int_comp_row(row_of_01,by=8):\n",
    "    ll = divide_row_by(row_of_01,by=by)\n",
    "    ss = []\n",
    "    for l in ll:\n",
    "        rl = run_length_row(l)\n",
    "        ss.append([v for k,v in rl])\n",
    "    return list(ss)    \n",
    "\n",
    "def divide_row_by(row,by):\n",
    "    ll = []\n",
    "    n = len(row)\n",
    "    m = n//by\n",
    "    #print(m)\n",
    "    for k in range(m):\n",
    "        ll.append(row[k*by:((k+1)*by)])\n",
    "    return list(ll)    \n",
    "\n",
    "def generateNotes(pitchlist,alphaPitch={\"pitch\":1,\"volume\":2,\"pause\":3},shuffle_notes=True):\n",
    "    from itertools import product\n",
    "    from music21 import pitch\n",
    "    #pitchlist = [p for p in list(range(60-1*octave*12,60+24-1*octave*12))]\n",
    "    #distmat = np.array(matrix([[np.sqrt(2*(1.0-kernPitch(x,y))) for x in pitchlist] for y in pitchlist]))\n",
    "    #permutation,distance = tspWithDistanceMatrix(distmat,exact=False)\n",
    "    #pitchlist = [pitchlist[permutation[k]] for k in range(len(pitchlist))]\n",
    "    print([pitch.Pitch(midi=int(p)) for p in pitchlist])\n",
    "    #durationlist = [n for n in durs]\n",
    "    #if len(durs)>2:\n",
    "    #    distmat = np.array(matrix([[np.sqrt(2*(1.0-kernDuration(x,y))) for x in durationlist] for y in durationlist]))\n",
    "    #    permutation,distance = tspWithDistanceMatrix(distmat)\n",
    "    #    durationlist = [durationlist[permutation[k]] for k in range(len(durationlist))]\n",
    "    #print(durationlist)\n",
    "    volumelist = vols = [(128//8)*(k+1) for k in range(8)] #[x*127 for x in [1.0/6.0,1.0/3.0,1.0/2.0,2.0/3.0 ]]\n",
    "    print(volumelist)\n",
    "    #distmat = np.array(matrix([[np.sqrt(2*(1.0-kernVolume(x,y))) for x in volumelist] for y in volumelist]))\n",
    "    #permutation,distance = tspWithDistanceMatrix(distmat)\n",
    "    #volumelist = [volumelist[permutation[k]] for k in range(len(volumelist))]\n",
    "    print(volumelist)\n",
    "    pauselist = [False,True]\n",
    "    ll = list(product(pauselist,volumelist,pitchlist))\n",
    "    if shuffle_notes:\n",
    "        shuffle(ll)\n",
    "    #distmat = np.array(matrix([[distKern(x,y,alphaPitch) for x in ll] for y in ll]))\n",
    "    #np.random.seed(43)\n",
    "    #permutation,distance = tspWithDistanceMatrix(distmat,exact=False)\n",
    "    #ll = [ll[permutation[k]] for k in range(len(ll))]\n",
    "    print(len(ll))\n",
    "    #print(ll)\n",
    "    pitches = [p[2] for p in ll]\n",
    "    #durations = [d[0] for d in ll]\n",
    "    volumes = [v[1] for v in ll]\n",
    "    isPauses = [p[0] for p in ll]\n",
    "    #print(pitches)\n",
    "    return pitchlist,volumelist,pauselist\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",
    "def get_knn_model(for_list,kern=kernPitch):\n",
    "    #notes = np.array([[x*1.0 for x in n] for n in notes])\n",
    "    from sklearn.neighbors import NearestNeighbors\n",
    "    np.random.seed(0)\n",
    "    #nbrs = NearestNeighbors( algorithm='ball_tree',metric=distKern(kern(zz0,alphaPitch=alphaPitch))).fit([[r] for r in range(len(zz0))])\n",
    "    M = matrix([[float(kern(a,b)) for a in for_list] for b in for_list],ring=RDF)\n",
    "    #print(M)\n",
    "    Ch = np.array(M.cholesky())\n",
    "    nbrs = NearestNeighbors().fit(Ch)\n",
    "    return nbrs, Ch\n",
    "\n",
    "def findBestMatches(nbrs,new_row,n_neighbors=3):\n",
    "    distances,indices = nbrs.kneighbors([np.array(new_row)],n_neighbors=n_neighbors)\n",
    "    dx = sorted(list(zip(distances[0],indices[0])))\n",
    "    #print(dx)\n",
    "    indi = [d[1] for d in dx]\n",
    "    #print(indi)\n",
    "    #print(distances)\n",
    "    #distances,indices = nbrs.query([np.array(new_row)],k=n_neighbors)\n",
    "    return indi\n",
    "\n",
    "\n",
    "from itertools import groupby\n",
    "\n",
    "def get_flattened_durations_from_graph(G,permutation_of_vertices,sorted_reverse,by=16,shuffled=True):\n",
    "    zz = sorted(list(zip(G.degree_sequence(),G.vertices())),reverse=sorted_reverse)\n",
    "    #print(sorted(zz))\n",
    "    if shuffled:\n",
    "        shuffle(zz)\n",
    "    A = G.adjacency_matrix(vertices = permutation_of_vertices)\n",
    "    ll = []\n",
    "    for row in A:\n",
    "        #print(row)\n",
    "        ll.extend([xx/by for xx in x] for x in (int_comp_row(row,by=by)))\n",
    "    ss = []\n",
    "    for l in ll:\n",
    "        ss.extend(l)\n",
    "    return ss    \n",
    "\n",
    "def get_bitstring_from_durations(durs,denom=None):\n",
    "    # convert reciprocal into rational number:\n",
    "    # take the denominator\n",
    "    # take the gcd of all denominators\n",
    "    qq = ([QQ(d) for d in durs])\n",
    "    dd = [q.denominator() for q in qq]\n",
    "    D = denom if not denom is None else lcm(dd)\n",
    "    nn = [q.numerator()*D/q.denominator() for q  in qq]\n",
    "    xx = []\n",
    "\n",
    "    print(D,nn)\n",
    "    y = 0\n",
    "    for n in nn:\n",
    "        xx.extend(n*[y])\n",
    "        y = 1*(not y)\n",
    "    return xx,D    \n",
    "    \n",
    " \n",
    "def get_durations_from_graph(G,sorted_reverse,by=16,shuffled=True):\n",
    "    zz = sorted(list(zip(G.degree_sequence(),G.vertices())),reverse=sorted_reverse)\n",
    "    #print(sorted(zz))\n",
    "    if shuffled:\n",
    "        shuffle(zz)\n",
    "    A = G.adjacency_matrix(vertices = [z[1] for z in zz])\n",
    "    ss = []\n",
    "    for row in A:\n",
    "        print(row)\n",
    "        xx = ((int_comp_row(row,by=by)))\n",
    "        ll = []\n",
    "        for x in xx: \n",
    "            ll.extend([xs/by for xs in x])\n",
    "        ss.append(ll)  \n",
    "    return ss    \n",
    "\n",
    "\n",
    "def interval_graph(intervals):\n",
    "    \"\"\"Generates an interval graph for a list of intervals given.\n",
    "\n",
    "    In graph theory, an interval graph is an undirected graph formed from a set\n",
    "    of closed intervals on the real line, with a vertex for each interval\n",
    "    and an edge between vertices whose intervals intersect.\n",
    "    It is the intersection graph of the intervals.\n",
    "\n",
    "    More information can be found at:\n",
    "    https://en.wikipedia.org/wiki/Interval_graph\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    intervals : a sequence of intervals, say (l, r) where l is the left end,\n",
    "    and r is the right end of the closed interval.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    G : networkx graph\n",
    "\n",
    "    Examples\n",
    "    --------\n",
    "    >>> intervals = [(-2, 3), [1, 4], (2, 3), (4, 6)]\n",
    "    >>> G = nx.interval_graph(intervals)\n",
    "    >>> sorted(G.edges)\n",
    "    [((-2, 3), (1, 4)), ((-2, 3), (2, 3)), ((1, 4), (2, 3)), ((1, 4), (4, 6))]\n",
    "\n",
    "    Raises\n",
    "    ------\n",
    "    :exc:`TypeError`\n",
    "        if `intervals` contains None or an element which is not\n",
    "        collections.abc.Sequence or not a length of 2.\n",
    "    :exc:`ValueError`\n",
    "        if `intervals` contains an interval such that min1 > max1\n",
    "        where min1,max1 = interval\n",
    "    \"\"\"\n",
    "    graph = nx.Graph()\n",
    "\n",
    "    tupled_intervals = [tuple(interval) for interval in intervals]\n",
    "    graph.add_nodes_from(tupled_intervals)\n",
    "\n",
    "    while tupled_intervals:\n",
    "        min1, max1,note1 = interval1 = tupled_intervals.pop()\n",
    "        graph.add_edge(interval1, interval1)\n",
    "        for interval2 in tupled_intervals:\n",
    "            min2, max2,note2 = interval2\n",
    "            if kernInt(interval1,interval2)>0:\n",
    "                graph.add_edge(interval1, interval2)\n",
    "    return graph\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def getCoordinatesOf(intList,kernel=kernInt,nDim=None):\n",
    "    M0 = matrix([[kernel(t1,t2) for t1 in intList] for t2 in intList],ring=RDF)\n",
    "    #print(M0)\n",
    "    from sklearn.decomposition import PCA\n",
    "    from sklearn.decomposition import KernelPCA\n",
    "    from sklearn.preprocessing import MinMaxScaler,StandardScaler\n",
    "    scaler = StandardScaler() #MinMaxScaler((0,1))\n",
    "    KPCA = KernelPCA(n_components=nDim,kernel='precomputed',eigen_solver='randomized')\n",
    "    \n",
    "    Ch0 = KPCA.fit_transform(np.array(M0))\n",
    "    #print(Ch0)\n",
    "    X0 = [x for x in 1.0*Ch0]    \n",
    "    \n",
    "    #print(X0)\n",
    "        \n",
    "    #X0 = scaler.fit_transform(X0)\n",
    "    \n",
    "    invPitchDict = dict(zip(intList,range(len(intList))))\n",
    "    return Ch0, invPitchDict\n",
    "\n",
    "def run_length_row(row_of_01):\n",
    "    from itertools import groupby\n",
    "    ar = row_of_01\n",
    "    return [(k, sum(1 for i in g)) for k,g in groupby(ar)]\n",
    "\n",
    "def pivoted_cholesky(M,max_rank=2):\n",
    "    import tensorflow_probability as tfp\n",
    "    import tensorflow as tf\n",
    "    import numpy as np\n",
    "    piv_chol = (tfp.math.pivoted_cholesky(M,max_rank = max_rank, diag_rtol=float(0.0001), name=None)).numpy()\n",
    "    return piv_chol\n",
    "\n",
    "def gco(il,kern,nDim=None):\n",
    "    M0 = matrix([[kern(t1,t2) for t1 in il] for t2 in il],ring=QQ)\n",
    "    from sklearn.preprocessing import MinMaxScaler,StandardScaler\n",
    "    scaler = StandardScaler() #MinMaxScaler((0,1))\n",
    "    from sklearn.decomposition import PCA\n",
    "    pca = PCA(n_components=nDim)\n",
    "    #KPCA = KernelPCA(n_components=nDim,kernel='precomputed',eigen_solver='randomized')\n",
    "    \n",
    "    Ch0 = pivoted_cholesky(M0,max_rank=nDim)\n",
    "    #print(Ch0)\n",
    "    X0 = [x for x in 1.0*Ch0]    \n",
    "    \n",
    "    #print(X0)\n",
    "        \n",
    "    X0 = pca.fit_transform(X0)\n",
    "    \n",
    "    invPitchDict = dict(zip(il,range(len(il))))\n",
    "    return X0, invPitchDict\n",
    "\n",
    "def get_knn_model(X):\n",
    "    #notes = np.array([[x*1.0 for x in n] for n in notes])\n",
    "    from sklearn.neighbors import NearestNeighbors\n",
    "    np.random.seed(0)\n",
    "    nbrs = NearestNeighbors( algorithm='ball_tree').fit(X)\n",
    "    return nbrs\n",
    "\n",
    "def findBestMatches(nbrs,new_row,n_neighbors=1):\n",
    "    distances,indices = nbrs.kneighbors([np.array(new_row)],n_neighbors=n_neighbors)\n",
    "    dx = sorted(list(zip(distances[0],indices[0])))\n",
    "    #print(dx)\n",
    "    indi = [d[1] for d in dx]\n",
    "    #print(indi)\n",
    "    #print(distances)\n",
    "    #distances,indices = nbrs.query([np.array(new_row)],k=n_neighbors)\n",
    "    return indi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9584cfb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "intervals = [(0,1/4),(1/4,2/4),(2/4,3/4),(3/4,4/4), # first voice\n",
    "             (0,1/2),(1/2,1/2+1/4),(3/4,1), # second voice\n",
    "             #(0,1/2),(1/2,1), # third voice\n",
    "             #(0,1/8),(1/8,1/2),(1/2,1/2+1/8),(1/2+1/8,1)\n",
    "            ]\n",
    "#G = nx.interval_graph(intervals)\n",
    "#sorted(G.edges)\n",
    "\n",
    "#nx.draw(G,with_labels=True)\n",
    "\n",
    "\n",
    "files = [\n",
    "    #\"./bsp.mid\",\n",
    "    #\"./bsp2.mid\",\n",
    "    #\"./melodie_floete.mid\",\n",
    "\"./una_mattina.mid\",\n",
    "\"./for_elise_by_beethoven.mid\",\n",
    "\"./For_You.mid\",\n",
    "#\"./Try_Again.mid\",\n",
    "#\"./Fly_Einaudi.mid\",\n",
    "\"./Piano_and_Cello.mid\",\n",
    "#\"./midi/Gymnopdie_No._1__Satie.mid\",\n",
    "#\"./midi/text_8_9.mid\",\n",
    "#\"./midi/experiment-Cyclic-group-of-order-4-as-a-permutation-groupgraph_48-8_4.mid\",\n",
    "    \"./midi/experiment_for_sine_waves.mid\",\n",
    "#    \"./midi/p-adic_numbers_and_partitions.mid\",\n",
    "#    \"./midi/algorithmic_composition_bridge_tower_limburg.mid\"\n",
    "]\n",
    "\n",
    "tempos = [\n",
    "    64,\n",
    "    70,\n",
    "    80,\n",
    "    80,\n",
    "    120\n",
    "]\n",
    "\n",
    "c = 0\n",
    "for file in files:\n",
    "    print(file)\n",
    "    scores = parse_file(file)\n",
    "#scores = parse_file(\"./for_elise_by_beethoven.mid\")\n",
    "#scores = parse_file(\"./For_You.mid\")\n",
    "#scores = parse_file(\"./Try_Again.mid\")\n",
    "#scores = parse_file(\"./Fly_Einaudi.mid\")\n",
    "#scores = parse_file(\"./Piano_and_Cello.mid\")\n",
    "#scores = parse_file(\"./midi/Gymnopdie_No._1__Satie.mid\")\n",
    "#scores = parse_file(\"./midi/text_8_9.mid\")\n",
    "    ints = []\n",
    "    notes = []\n",
    "    dintpitch = dict([])\n",
    "    iinds = []\n",
    "    notesSet = set([])\n",
    "    nDim = None\n",
    "    \n",
    "    for sc in scores:\n",
    "        t = 0\n",
    "        notes = []\n",
    "        notesSet = set([])\n",
    "        ints = []\n",
    "        for note in sc[0:]:\n",
    "            pitch,duration,volume,isPause = note\n",
    "            notes.append(tuple(note))\n",
    "            notesSet.add(tuple(note))\n",
    "            #if not isPause:\n",
    "            ints.append((t,t+duration,tuple(note)))\n",
    "                #dintpitch[(t,t+duration)] = note\n",
    "            t += duration\n",
    "            \n",
    "#print(ints)    \n",
    "        lns = list(notesSet)\n",
    "        print(len(ints))\n",
    "        #nDim = min(len(ints),100)\n",
    "        XNotes , xnotes = getCoordinatesOf(ints,kernel=kernInt,nDim=nDim)\n",
    "        #print(XNotes)\n",
    "  \n",
    "        import scipy, numpy as np\n",
    "\n",
    "        nbrsNotes = get_knn_model(XNotes)\n",
    "    \n",
    "        #XXPitch = np.array([XPitch[xpitch[note[0]]] for note in notes])\n",
    "        #XXDuration = np.array([XDuration[xduration[note[1]]] for note in notes])\n",
    "        #XXVolume = np.array([XVolume[xvolume[note[2]]] for note in notes])\n",
    "        #XXNotes = np.array([XNotes[xnotes[note]] for note in notes])\n",
    "        #XX = np.array([(XPitch[xpitch[note[0]]][0],XDuration[xduration[note[1]]][0],XVolume[xvolume[note[2]]][0]) for note in notes])\n",
    "    #pg = plot(matrix_plot(matrix(np.round(X,2))))\n",
    "\n",
    "        #XBeziers = [XXPitch,XXDuration,XXVolume]\n",
    "        XBeziers = [XNotes]\n",
    "        nbrs = [nbrsNotes]\n",
    "        dicts = [xnotes]\n",
    "    \n",
    "        new_notes = []\n",
    "        div = 30*len(ints)\n",
    "        #for s in [i/div*1.0 for i in range(0,div+1)]:\n",
    "        new_sth = []\n",
    "        for i in range(len(XBeziers)):\n",
    "            import numpy as np , bezier\n",
    "            XXBezier = XBeziers[i]\n",
    "            nodes = XXBezier.T\n",
    "            curve = bezier.Curve(nodes, degree=XXBezier.shape[0]-1)\n",
    "            ss = [np.sin((j/div+0.5)*2*np.pi)/2.0+0.5 for j in range(0,div+1)]\n",
    "            ss = [np.cos((j/div)*2*np.pi)/2.0+0.5 for j in range(0,div+1)]\n",
    "            #ss = [(j/div*1.0)**3 for j in range(0,div+1)]\n",
    "            #ss = [(np.exp(j)/np.exp(div)*1.0)**2 for j in range(0,div+1)]\n",
    "            print(ss)\n",
    "    #plt.show(curve.plot(num_pts=X.shape[0]))\n",
    "            for s in ss:\n",
    "                new_row = curve.evaluate(s)\n",
    "                #print(s,new_row)\n",
    "            #print(new_row.T[0])\n",
    "                indi = findBestMatches(nbrs[i],new_row.T[0])\n",
    "                #print(indi)\n",
    "                #print(ll[i])\n",
    "                new_notes.append(ints[int(indi[0])][2])\n",
    "            \n",
    "        #mynotes = list(zip(*new_notes))\n",
    "        #print(mynotes)\n",
    "        mmnotes = [r[0] for r in run_length_row(new_notes)]\n",
    "        #print(mmnotes)\n",
    "        print(new_notes)\n",
    "        iinds.append(mmnotes)\n",
    "    #pg.save(\"matrix_plot_\"+file+\".png\")\n",
    "    #plt.plot(pg)\n",
    "    #plt.show(pg)\n",
    "    #print(np.round(X,2))\n",
    "    #A = ((np.round(np.dot(X,X.T),12)>0)*1)\n",
    "    #plot(Graph(nx.Graph(A),loops=True),**{\"vertex_labels\":False,\"layout\":\"spring\",\"vertex_size\":0.2,\"loop_size\":0.01}).show()\n",
    "    #print((np.round(np.dot(X,X.T),2)))\n",
    "\n",
    "#    print(invIntDict)\n",
    "\n",
    "    #H = interval_graph(ints)\n",
    "    #H = Graph(G,loops=True)\n",
    "    #sorted(G.edges)\n",
    "    \n",
    "    #m = matrix([[kernInt(i1,i2) for i1 in H.nodes] for i2 in H.nodes],ring=RDF)\n",
    "    #print(m.is_positive_definite())\n",
    "    #print(m.is_symmetric())\n",
    "    #print(sum(sum(m))/len(H.edges()))\n",
    "    \n",
    "    #print(sorted(m.eigenvalues()))\n",
    "\n",
    "    #nx.draw_spectral(H,with_labels=False)\n",
    "\n",
    "            #def ff(s,i=0):\n",
    "            #    x = curve.evaluate(s)\n",
    "            #    return x[i,0]\n",
    "\n",
    "            #parametric_plot3d((lambda s: ff(s,i=0),lambda s: ff(s,i=1),lambda s: ff(s,i=2)), (0,1),plot_points=1000).show()\n",
    "    \n",
    "    #G = Graph(H,loops=True)\n",
    "#print(sorted(scipy.linalg.eigh(matrix(G.laplacian_matrix(),ring=RDF).numpy())[0]))\n",
    "    #ccs = G.connected_components_sizes()\n",
    "    #giant = G.connected_components_subgraphs()[0]\n",
    "    #p = (len(G.vertices()),\n",
    "    #     len(G.edges()),\n",
    "    #     mean(G.degree_sequence())*1.0,\n",
    "    #     ccs[0]/sum(ccs)*1.0,\n",
    "    #     G.average_distance()*1.0,\n",
    "    #     #G.clustering_average()*1.0,\n",
    "    #     giant.average_distance()*1.0,\n",
    "    #     len([c for c in ccs if c>1])\n",
    "    #    )\n",
    "    #print(p)\n",
    "    #plot(G,**{\"vertex_labels\":False,\"layout\":\"spring\",\"vertex_size\":0.2,\"loop_size\":0.01}).show()\n",
    "    fn = file+\"-bezier.mid\"\n",
    "    writePitches(fn=fn,inds=iinds,tempo=tempos[c],instrument=len(iinds)*[0],add21=False,start_at= len(iinds)*[0],durationsInQuarterNotes=False)\n",
    "    c+=1\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ca3ac95",
   "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": 5
}