Hierarchical clustering 1: Dissimilarity Matrix

//---a helper function, similar to imagetransform for 1/2 dim waves
threadsafe Function/wave  ImageTransformRowCol(variable row,variable col,wave inputwave)
    make/o/n=1 ExtractedWave
    if(numpnts(inputwave)>1)
        if(numtype(row)==0)
            if((dimsize(inputwave,0)<=1))
                duplicate/o inputwave,ExtractedWave
            elseif((dimsize(inputwave,0)>1)&&(dimsize(inputwave,1)>0))
                imagetransform/g=(row) getrow inputwave
                wave W_ExtractedRow
                duplicate/o W_ExtractedRow,ExtractedWave
            endif
        else
            if((dimsize(inputwave,1)<=1))
                duplicate/o inputwave,ExtractedWave
            elseif((dimsize(inputwave,1)>1)&&(dimsize(inputwave,1)>col))
                imagetransform/g=(col) getcol inputwave
                wave W_ExtractedCol
                duplicate/o W_ExtractedCol,ExtractedWave
            endif
        endif
    endif
    return ExtractedWave
end
//---actual clustering
//---find the difference between two waves
Threadsafe function ComputeDissimilarity2waves(wave Cluster_Waves,variable  pos1,variable  pos2)
    duplicate/o ImageTransformRowCol(nan,pos1,Cluster_Waves),waveI
    duplicate/o ImageTransformRowCol(nan,pos2,Cluster_Waves),waveJ
    waveI=(waveI-waveJ)^2  //--- squared difference between the two waves
    return sqrt(sum(waveI))/numpnts(waveI)
end

//---creates the dissimilarity matrix based on the selected experiments
//---we assume that the wave Cluster_Waves is a 2D wave in which the columns represent different samples (for example, a recorded time-trace).

Function MakeDisSimilarityMatrix(wave Cluster_Waves)
    variable i,j,reg
    make/o/n=(dimsize(Cluster_Waves,1),dimsize(Cluster_Waves,1)) DisSimilarityMatrix=nan
    Variable n_,nthreads= ThreadProcessorCount
    Variable threadGroupID= ThreadGroupCreate(nthreads)    
    for(i=0;i<dimsize(DisSimilarityMatrix,0)-1;i++)
        j=i+1
        do
            for(n_=0;n_<nthreads;n_++)
                if(n_+j<dimsize(DisSimilarityMatrix,0))
                    ThreadStart threadGroupID,n_,ComputeDissimilarity2waves(Cluster_Waves,i,n_+j)//---populate the similarity matrix
                endif
            endfor
            do
                Variable threadGroupStatus = ThreadGroupWait(threadGroupID,1000)
            while( threadGroupStatus != 0 )//---release threads        
            for(n_=0;n_<nthreads;n_++)
                if(n_+j<dimsize(DisSimilarityMatrix,0))
                    DisSimilarityMatrix[i][n_+j]=ThreadReturnValue(threadGroupID,n_)
                    DisSimilarityMatrix[n_+j][i]=DisSimilarityMatrix[i][n_+j]
                endif
            endfor                         
            j+=n_          
        while(j<dimsize(DisSimilarityMatrix,0))
    endfor
    Variable dummy= ThreadGroupRelease(-2)  //---release threads
end