Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • RMTGeneNET, c program is not running

    /* v2: Replaces the GSL eigenvalue solver with the MAGMA solver */
    /* v2-1: Branching from v2, attempting to use the ACML lapack
    * eigenvector solver instead of MAGMA (or perhaps as an
    * intermediate step. */
    /*
    * v2-2: Uses MKL rather than ACML for eigen solve */
    /*
    * v2-2.1: Removed the sort call and eigen print inside of chiSquareTest */
    /*
    * v2-3: Uses a symmettric matrix eigen solver. Use with RMM.v1-3 */

    #include <stdlib.h>
    #include <stdio.h>
    #include <string.h>
    #include <gsl/gsl_interp.h>
    #include <gsl/gsl_spline.h>
    #include <gsl/gsl_math.h>
    #include <gsl/gsl_eigen.h>
    #include <gsl/gsl_matrix.h>
    #include <mkl.h>

    void swapD(double* l, int idx1, int idx2){
    double temp = l[idx1];
    l[idx1] = l[idx2];
    l[idx2] = temp;
    return;
    }

    void swapF(float* l, int idx1, int idx2){
    float temp = l[idx1];
    l[idx1] = l[idx2];
    l[idx2] = temp;
    return;
    }

    void quickSortD(double* l, int size){
    if(size<=1) return;
    int pivIdx = (int) size/1.618;//golden ratio
    double pivot = l[pivIdx];
    swapD(l, pivIdx, size-1);
    int leftPlace = 0;
    int i;
    for(i=0;i<size-1;i++){
    if(l[i]<pivot){
    swapD(l, i, leftPlace);
    leftPlace++;
    }
    }
    swapD(l, size-1, leftPlace);
    quickSortD(l,leftPlace);
    quickSortD(&l[leftPlace+1], size-leftPlace-1);
    return;
    }

    void quickSortF(float* l, int size){
    if(size<=1) return;
    int pivIdx = (int) size/1.618;//golden ratio
    float pivot = l[pivIdx];
    swapF(l, pivIdx, size-1);
    int leftPlace = 0;
    int i;
    for(i=0;i<size-1;i++){
    if(l[i]<pivot){
    swapF(l, i, leftPlace);
    leftPlace++;
    }
    }
    swapF(l, size-1, leftPlace);
    quickSortF(l,leftPlace);
    quickSortF(&l[leftPlace+1], size-leftPlace-1);
    return;
    }

    float* calculateEigen(float* mat, int size){
    float* W = (float*) malloc(sizeof(float)*size);
    float* work = (float*) malloc(sizeof(float)*5*size);
    MKL_INT rc;
    char jobz = 'N';//don't compute eigenvectors
    char uplo = 'U';//upper is stored (doesn't matter, both are stored)
    MKL_INT lwork = 5*size;
    //ssyev( char *jobz, char *uplo, MKL_INT *n, float *a, MKL_INT *lda, float *w, float *work, MKL_INT *lwork, MKL_INT *info );
    ssyev(&jobz, &uplo , &size, mat, &size, W, work, &lwork ,&rc);
    if(rc!=0){
    printf("\nThe SSYEV function encountered a problem while solving, error code %d. The RandomMatrixModeling will continue without interuption, however it may be important to note this error.\n", rc);
    }
    free(work);
    return W;
    }

    //returned array will always be sorted and of length size-1
    double* unfolding(float* e, int size, int m){
    int count=1, i,j=0;//count equals 1 initially because of 2 lines following loop which propogates the arrays
    for(i=0; i<size-m; i+=m) count++;
    double* oX = (double*) malloc(sizeof(double)*count);
    double* oY = (double*) malloc(sizeof(double)*count);
    for(i=0; i<size-m; i+=m){
    oX[j]=e[i];
    oY[j]= (i+1.0)/(double)size;
    j++;
    }
    oX[count-1] = e[size-1];
    oY[count-1] = 1;

    for(i=1;i<count;i++){
    if(!(oX[i-1]<oX[i])){
    printf("\nat postion %d a problem exists\n", i);
    printf("oX[i-1]=%f whilst oX[i]=%f\n",oX[i-1],oX[i]);
    }
    }
    double* yy = (double*) malloc(sizeof(double)*(size));

    gsl_interp_accel *acc = gsl_interp_accel_alloc();
    gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline, count);//see gsl docs, chapter 27: cspline is a natural spline
    gsl_spline_init(spline, oX, oY, count);

    for(i=0;i<(size-2);i++){
    yy[i+1] = gsl_spline_eval(spline, e[i+1], acc);
    }
    gsl_spline_free(spline);
    gsl_interp_accel_free(acc);
    yy[0] = 0.0;
    yy[size-1] = 1.0;
    for(i=0;i<size-1;i++){
    yy[i] = size*(yy[i+1]-yy[i]);
    }
    quickSortD(yy, size-1);
    free(oX);
    free(oY);
    return yy;
    }

    float* degenerate(float* eigens, int size, int* newSize){
    int i, j=0, count=1;//because one flag is set before the loop
    for(i=0;i<size;i++){
    if(fabs(eigens[i]) < 0.000001){
    eigens[i] = 0.0;
    }
    }
    int* flags = (int*) malloc(sizeof(int)*size);
    memset(flags, 0, size*sizeof(int));
    float temp = eigens[0];
    flags[0]=1;
    for(i=1;i<size; i++){
    if(fabs(eigens[i]-temp) > 0.000001){
    count++;
    flags[i] = 1;
    temp = eigens[i];
    }
    }
    float* remDups = (float*) malloc(sizeof(float)*count);//remDups means "removed duplicates"
    for(i=0;i<size;i++){
    if(flags[i]==1){
    remDups[j] = eigens[i];
    j++;
    }
    }
    free(flags);
    *newSize = count;

    return remDups;
    }


    double chiSquareTestUnfoldingNNSDWithPoisson4(float* eigens, int size, double bin, int pace){
    int newSize;
    float* newE;
    double* edif;
    newE = degenerate(eigens, size, &newSize);
    size = newSize;

    edif = unfolding(newE, size, pace);
    free(newE);
    size = size-1; //see note above unfolding function, will return an array of size-1
    int n = (int) (3.0/bin) + 1;
    double obj, expect, chi = 0;
    int i, j, count;
    for(i=0;i<n;i++){
    count = 0;
    for(j=0;j<size;j++){
    if(edif[j]>i*bin && edif[j] < (i+1)*bin) count++;
    }
    obj = (double) count;
    expect = (exp(-1*i*bin)-exp(-1*(i+1)*bin))*size;
    chi += (obj -expect)*(obj-expect)/expect;
    }
    free(edif);
    return chi;
    }

    //calls same name, 4 args instead of 5
    double chiSquareTestUnfoldingNNSDWithPoisson(float* eigens, int size, double bin, int minPace, int maxPace){
    double chiTest =0;
    int i = 0;
    int m;

    i=0;
    for(m = minPace; m<maxPace; m++){
    chiTest+=chiSquareTestUnfoldingNNSDWithPoisson4(eigens, size, bin, m);
    i++;
    }

    return chiTest/i;
    }

    this programm is not running,anybody help me please?

Latest Articles

Collapse

  • seqadmin
    The Impact of AI in Genomic Medicine
    by seqadmin



    Artificial intelligence (AI) has evolved from a futuristic vision to a mainstream technology, highlighted by the introduction of tools like OpenAI's ChatGPT and Google's Gemini. In recent years, AI has become increasingly integrated into the field of genomics. This integration has enabled new scientific discoveries while simultaneously raising important ethical questions1. Interviews with two researchers at the center of this intersection provide insightful perspectives into...
    02-26-2024, 02:07 PM
  • seqadmin
    Multiomics Techniques Advancing Disease Research
    by seqadmin


    New and advanced multiomics tools and technologies have opened new avenues of research and markedly enhanced various disciplines such as disease research and precision medicine1. The practice of merging diverse data from various ‘omes increasingly provides a more holistic understanding of biological systems. As Maddison Masaeli, Co-Founder and CEO at Deepcell, aptly noted, “You can't explain biology in its complex form with one modality.”

    A major leap in the field has
    ...
    02-08-2024, 06:33 AM

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, 02-23-2024, 04:11 PM
0 responses
57 views
0 likes
Last Post seqadmin  
Started by seqadmin, 02-21-2024, 08:52 AM
0 responses
67 views
0 likes
Last Post seqadmin  
Started by seqadmin, 02-20-2024, 08:57 AM
0 responses
56 views
0 likes
Last Post seqadmin  
Started by seqadmin, 02-14-2024, 09:19 AM
0 responses
65 views
0 likes
Last Post seqadmin  
Working...
X