Violin plots

#pragma rtGlobals=3     // Use modern global access method and strict wave access.
#include <Percentile and Box Plot>
 
//Before executing, you need to have plotted the boxplots using Percentile and Box Plot
//Plot the box plots into a target window using all waves from a source window
//This is a limitation of the functions, they rely on two graph windows (source and target) for making the violins.
//A width of 0.1 is good for the box plots and the boxplots look better with whiskers going to 0 and 100 rather than 10 and 90 (Igor default)
 
Menu "Macros"
    "Violin Plot...",  ViolinPlot()
    "Scrub Violins", ClearViolins()
End
 
Function ViolinPlot ()
    String sourceWindow // window with the waves to make boxplot and violins
    String targetWindow // window with the boxplots to add violins
    Variable trimValue=1    // violins can be trimmed (default) or plotted with long tails
    Variable thickValue=0.1 // this is the width of the boxplot
    
    Prompt sourceWindow, "What is the source graph?", popup, WinList("*", ";", "WIN:1")
    Prompt targetWindow, "What is the box plot graph?", popup, WinList("*", ";", "WIN:1")
    Prompt  trimValue, "yes=1, no=0"
    Prompt thickValue, "auto=0.1"
    DoPrompt "Pick graphs", sourceWindow, targetWindow, trimValue, thickValue
    
    String/G gSrcWin = sourceWindow //set global strings and variables
    String/G gTgtWin = targetWindow
    Variable/G gtrim = trimValue
    Variable/G gthick = thickValue
    
    DoWindow /F $sourceWindow
    String wlist = Wavelist("*",";","WIN:")
    String name
    Variable nWaves = ItemsInList(wlist)
    
    Variable i
    
    for(i = 0; i < nWaves; i += 1)
        name = StringFromList(i,wList)  //picks waves from source window and sends them to be plotted as violins in target window
        Wave w0 = $name
        Violin(w0,i)
    endfor
    
    DoWindow /F $targetWindow
    SetAxis/A/N=1/E=1 left
    SetAxis bottom -0.5,(nWaves-0.5)
End
 
///     @param  w   1D wave for plotting as violin
///     @param  xpos    variable to describe where violin will go on x-axis
Function Violin(w,xpos)
    Wave w
    Variable xpos
    
    SVAR srcWin = gSrcWin //set local ref for global strings
    SVAR tgtWin = gTgtWin
    NVAR trim = gtrim
    NVAR thick = gThick
    
    Variable nPoints = numpnts(w)
    Variable xMin, xMax
    
    if(trim == 1)
        xMin = wavemin(w)
        xMax = wavemax(w)                       // trimmed violins
    else
        Variable bw = ((4 * sqrt(variance(w))^5) / (3 * numpnts(w)))^0.2    // this is Silverman's rule-of-thumb
        xMin = wavemin(w) - (2*bw)
        xMax = wavemax(w) + (2*bw)              // not trimmed, limit kde to 2 bandwidths beyond min/max
    endif
 
    Variable dX = min((xMax - xMin) / 100,1) // delta x for StatsKDE
    StatsKDE/Q/BWM=2/DEST=kdeWave/KT=1/S={xMin,dX,xMax} w
    WAVE/Z w1 = kdeWave
    Variable auc = sum(w1)
    w1 /= auc * (thick / 2) // omit *thick for a violin=2 i.e. two half-violins that =1
    
    String wnL = NameofWave(w) + "_kdeL"
    Duplicate/D/O w1 $wnL
    String wnR = NameofWave(w) + "_kdeR"
    Duplicate/D/O w1 $wnR
    String wnLx = NameofWave(w) + "_kdeLx"
    Duplicate/D/O w1 $wnLx
    String wnRx = NameofWave(w) + "_kdeRx"
    Duplicate/D/O w1 $wnRx
    String vwny = NameofWave(w) + "_kdevy"
    String vwnx = NameofWave(w) + "_kdevx"
    KillWaves/Z w1
    Wave wnL1 = $wnL
    Wave wnR1 = $wnR
    Wave wnLx1 = $wnLx
    Wave wnRx1 = $wnRx
    wnLx1 = x                   //get x scale
    wnRx1 = x
    
    if (trim == 1)
        InsertPoints 0,1, wnLx1 //add point, set to same as first point to complete the bottom of a trimmed violin
        InsertPoints 0,1, wnRx1
        wnLx1[0] = wnLx1[1]
        wnRx1[0] = wnRx1[1]
        InsertPoints 0,1, wnL1
        InsertPoints 0,1, wnR1
        wnL1[0] = 0
        wnR1[0] = 0
    endif
    
    wnR1 *= -1                                      //make the other side of violin
    Reverse wnR1, wnRx1                             //reverse the righthand side of violin
    Concatenate/O/NP {wnL1,wnR1}, $vwny     // form the wave for y position of outline (will be x when vertical)
    Concatenate/O/NP {wnLx1,wnRx1}, $vwnx       // form the wave for x position of outline (will be y when vertical)
    wave vwny1 = $vwny
    wave vwnx1 = $vwnx
    
    DoWindow/F $tgtWin                                                      //bring target window to the front
    SetDrawLayer ProgBack                                                       // violins are drawn in ProgBack
    SetDrawEnv ycoord=left, xcoord=bottom, fillfgc=(65535,0,0),fillpat=3        // violins are red
    DrawPoly xpos+vwny1[0],vwnx1[0],1,1, vwny1,vwnx1                    // draws violin
    KillWaves/Z wnL1,wnR1,wnLx1,wnRx1                                       // clean-up
    SetDrawLayer UserFront                                                  // set layer back
End
 
//Handy function to scrub the violins away
Function ClearViolins()
    SVAR tgtWin = gTgtWin
    DoWindow /F $tgtWin
    SetDrawLayer /K ProgBack // kills ProgBack
    SetDrawLayer UserFront
End