--- /dev/null
+CC = gcc
+CFLAGS = -g -std=gnu99 -Wno-implicit-function-declaration
+LDFLAGS = -Lobj/ -lvalse_core
+LIB = valse_core.so
+LIB_SRC = $(wildcard ../sources/*.c)
+LIB_OBJ = $(LIB_SRC:.c=.o)
+
+all: $(LIB) test.EMGLLF test.EMGrank test.constructionModelesLassoMLE test.EMGrank test.constructionModelesLassoRank test.selectionTotale
+
+$(LIB): $(LIB_OBJ)
+ $(CC) -o $@ $^
+
+test.EMGLLF: test.EMGLLF.o
+ $(CC) -o $@ $^ $(LDFLAGS)
+
+test.constructionModelesLassoMLE: test.constructionModelesLassoMLE.o
+ $(CC) -o $@ $^ $(LDFLAGS)
+
+test.EMGrank: test.EMGrank.o
+ $(CC) -o $@ $^ $(LDFLAGS)
+
+test.constructionModelesLassoRank: test.constructionModelesLassoRank.o
+ $(CC) -o $@ $^ $(LDFLAGS)
+
+test.selectionTotale: test.selectionTotale.o
+ $(CC) -o $@ $^ $(LDFLAGS)
+
+%.o: %.c
+ $(CC) -o $@ -c $< $(CFLAGS)
+
+clean:
+ rm -f *.o ../sources/*.o
+
+cclean: clean
+ rm -f $(LIB) $(TEST)
+
+.PHONY: all clean cclean
--- /dev/null
+1) transformer les .m en .R
+2) sauvegarder les résultats + entrée/sorties à partir du code R uniquement
+3) dans test.Truc.c, prendre en entrée les paramètres de data/ juste sauvegardés,
+ puis comparer les sorties aux sorties enregistrées
--- /dev/null
+function[phiInit,rhoInit,piInit,gamInit] = basicInitParameters(n,p,m,k)
+
+ phiInit = zeros(p,m,k);
+
+ piInit = (1.0/k) * ones(1,k);
+
+ rhoInit = zeros(m,m,k);
+ for r=1:k
+ rhoInit(:,:,r) = eye(m,m);
+ end
+
+ gamInit = 0.1 * ones(n,k);
+ R = random('unid',k,n,1);
+ for i=1:n
+ gamInit(i,R(i)) = 0.9;
+ end
+ gamInit = gamInit / (sum(gamInit(1,:)));
+
+end
--- /dev/null
+function[]=checkOutput(varName, matrix, refMatrix, tol)
+
+ fprintf('Checking %s\n',varName);
+ maxError = max(max(max(max(abs(matrix - refMatrix)))));
+ if maxError >= tol
+ fprintf(' Inaccuracy: max(abs(error)) = %g >= %g\n',maxError,tol);
+ else
+ fprintf(' OK\n');
+ end
+
+end
--- /dev/null
+%covariance matrix for tests on synthetic data: A(i,j) = a^|i-j|
+function[A] = covariance(p,a)
+
+ A = a*ones(p,p);
+ for i=1:p
+ A(i,:) = A(i,:) .^ abs(i-(1:p));
+ end
+
+end
--- /dev/null
+%X is generated following a gaussian mixture \sum pi_r N(meanX_k, covX_r)
+%Y is generated then, with Y_i ~ \sum pi_r N(Beta_r.X_i, covY_r)
+function[X,Y,Z] = generateIO(meanX, covX, covY, pi, beta, n)
+
+ [p, ~, k] = size(covX);
+ [m, ~, ~] = size(covY);
+ if exist('octave_config_info')
+ %Octave statistics package doesn't have gmdistribution()
+ X = zeros(n, p);
+ Z = zeros(n);
+ cs = cumsum(pi);
+ for i=1:n
+ %TODO: vectorize ? http://stackoverflow.com/questions/2977497/weighted-random-numbers-in-matlab
+ tmpRand01 = rand();
+ [~,Z(i)] = min(cs - tmpRand01 >= 0);
+ X(i,:) = mvnrnd(meanX(Z(i),:), covX(:,:,Z(i)), 1);
+ end
+ else
+ gmDistX = gmdistribution(meanX, covX, pi);
+ [X, Z] = random(gmDistX, n);
+ end
+
+ Y = zeros(n, m);
+ BX = zeros(n,m,k);
+ for i=1:n
+ for r=1:k
+ %compute beta_r . X_i
+ BXir = zeros(1, m);
+ for mm=1:m
+ BXir(mm) = dot(X(i,:), beta(:,mm,r));
+ end
+ %add pi(r) * N(beta_r . X_i, covY) to Y_i
+ Y(i,:) = Y(i,:) + pi(r) * mvnrnd(BXir, covY(:,:,r), 1);
+ end
+ end
+
+end
--- /dev/null
+%call generateIO with default parameters (random means, covariances = identity, equirepartition)
+function[X,Y,Z] = generateIOdefault(n, p, m, k)
+
+ rangeX = 100;
+ meanX = rangeX * (1 - 2*rand(k, p));
+ covX = zeros(p,p,k);
+ covY = zeros(m,m,k);
+ for r=1:k
+ covX(:,:,r) = eye(p);
+ covY(:,:,r) = eye(m);
+ end
+ pi = (1/k) * ones(1,k);
+
+ %initialize beta to a random number of non-zero random value
+ beta = zeros(p,m,k);
+ for j=1:p
+ nonZeroCount = ceil(m*rand(1));
+ beta(j,1:nonZeroCount,:) = rand(nonZeroCount, k);
+ end
+
+ [X,Y,Z] = generateIO(meanX, covX, covY, pi, beta, n);
+
+end
--- /dev/null
+function[] = generateRunSaveTest_EMGLLF(n, p, m, k, mini, maxi, gamma, lambda, varargin)
+
+ %set defaults for optional inputs
+ optargs = {200 15 10 3 5 10 1.0 0.5};
+ %replace defaults by user parameters
+ optargs(1:length(varargin)) = varargin;
+ [n, p, m, k, mini, maxi, gamma, lambda] = optargs{:};
+ tau = 1e-6;
+ mini = int64(mini);
+ maxi = int64(maxi);
+
+ %Generate phiInit,piInit,...
+ [phiInit,rhoInit,piInit,gamInit] = basicInitParameters(n, p, m, k);
+
+ %Generate X and Y
+ [X, Y, ~] = generateIOdefault(n, p, m, k);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ testFolder = 'data/';
+ mkdir(testFolder);
+ delimiter = '\n';
+
+ %save inputs
+ dlmwrite(strcat(testFolder,'phiInit'), reshape(phiInit,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'rhoInit'), reshape(rhoInit,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'piInit'), piInit, delimiter);
+ dlmwrite(strcat(testFolder,'gamInit'), reshape(gamInit,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'mini'), mini, delimiter);
+ dlmwrite(strcat(testFolder,'maxi'), maxi, delimiter);
+ dlmwrite(strcat(testFolder,'gamma'), gamma, delimiter);
+ dlmwrite(strcat(testFolder,'lambda'), lambda, delimiter);
+ dlmwrite(strcat(testFolder,'X'), reshape(X,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'Y'), reshape(Y,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'tau'), tau, delimiter);
+ dlmwrite(strcat(testFolder,'dimensions'), [n,p,m,k], delimiter);
+
+ [phi,rho,pi,LLF,S] = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau);
+
+ %save output
+ dlmwrite(strcat(testFolder,'phi'), reshape(phi,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'rho'), reshape(rho,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'pi'), pi, delimiter);
+ dlmwrite(strcat(testFolder,'LLF'), LLF, delimiter);
+ dlmwrite(strcat(testFolder,'S'), reshape(S,1,[]), delimiter);
+
+end
--- /dev/null
+function[] = generateRunSaveTest_constructionModelesLassoMLE(n, p, m, k, mini, maxi, gamma, glambda, varargin)
+
+ %set defaults for optional inputs
+ optargs = {200 15 10 3 5 10 1.0 [0.0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,0.85,0.99]};
+ %replace defaults by user parameters
+ optargs(1:length(varargin)) = varargin;
+ [n, p, m, k, mini, maxi, gamma, glambda] = optargs{:};
+ tau = 1e-6;
+ seuil = 1e-15;
+ mini = int64(mini);
+ maxi = int64(maxi);
+ L = length(glambda);
+
+ %Generate phiInit,piInit,...
+ [phiInit,rhoInit,piInit,gamInit] = basicInitParameters(n, p, m, k);
+
+ %Generate X and Y
+ [X, Y, ~] = generateIOdefault(n, p, m, k);
+
+ A2 = zeros(p,m+1,L,'int64');
+ for i=1:L
+ A2(:,1,i) = 1:p;
+ A2(1:5,2,i) = 1:5;
+ end
+ A1 = zeros(p,m+1,L,'int64');
+ for i=1:L
+ A1(:,1,i) = 1:p;
+ A1(1:5,2,i) = 1:5;
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ testFolder = 'data/';
+ mkdir(testFolder);
+ delimiter = '\n';
+
+ %save inputs
+ dlmwrite(strcat(testFolder,'phiInit'), reshape(phiInit,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'rhoInit'), reshape(rhoInit,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'piInit'), piInit, delimiter);
+ dlmwrite(strcat(testFolder,'gamInit'), reshape(gamInit,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'mini'), mini, delimiter);
+ dlmwrite(strcat(testFolder,'maxi'), maxi, delimiter);
+ dlmwrite(strcat(testFolder,'gamma'), gamma, delimiter);
+ dlmwrite(strcat(testFolder,'glambda'), glambda, delimiter);
+ dlmwrite(strcat(testFolder,'X'), reshape(X,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'Y'), reshape(Y,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'seuil'), seuil, delimiter);
+ dlmwrite(strcat(testFolder,'tau'), tau, delimiter);
+ dlmwrite(strcat(testFolder,'A1'), reshape(A1,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'A2'), reshape(A2,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'dimensions'), [n,p,m,k,L], delimiter);
+
+ [phi,rho,pi,lvraisemblance] = constructionModelesLassoMLE(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2);
+
+ %save output
+ dlmwrite(strcat(testFolder,'phi'), reshape(phi,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'rho'), reshape(rho,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'pi'), reshape(pi,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'lvraisemblance'), reshape(lvraisemblance,1,[]), delimiter);
+
+end
--- /dev/null
+function[] = generateRunSaveTest_selectiontotale(n, p, m, k, mini, maxi, gamma, glambda, varargin)
+
+ %set defaults for optional inputs
+ optargs = {200 15 10 3 5 10 1.0 [0.0,0.01,0.02,0.03,0.05,0.1,0.2,0.3,0.5,0.7,0.85,0.99]};
+ %replace defaults by user parameters
+ optargs(1:length(varargin)) = varargin;
+ [n, p, m, k, mini, maxi, gamma, glambda] = optargs{:};
+ tau = 1e-6;
+ seuil = 1e-15;
+ mini = int64(mini);
+ maxi = int64(maxi);
+ L = length(glambda);
+
+ %Generate phiInit,piInit,...
+ [phiInit,rhoInit,piInit,gamInit] = basicInitParameters(n, p, m, k);
+
+ %Generate X and Y
+ [X, Y, ~] = generateIOdefault(n, p, m, k);
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ testFolder = 'data/';
+ mkdir(testFolder);
+ delimiter = '\n';
+
+ %save inputs
+ dlmwrite(strcat(testFolder,'phiInit'), reshape(phiInit,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'rhoInit'), reshape(rhoInit,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'piInit'), piInit, delimiter);
+ dlmwrite(strcat(testFolder,'gamInit'), reshape(gamInit,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'mini'), mini, delimiter);
+ dlmwrite(strcat(testFolder,'maxi'), maxi, delimiter);
+ dlmwrite(strcat(testFolder,'gamma'), gamma, delimiter);
+ dlmwrite(strcat(testFolder,'glambda'), glambda, delimiter);
+ dlmwrite(strcat(testFolder,'X'), reshape(X,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'Y'), reshape(Y,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'seuil'), seuil, delimiter);
+ dlmwrite(strcat(testFolder,'tau'), tau, delimiter);
+ dlmwrite(strcat(testFolder,'dimensions'), [n,p,m,k,L], delimiter);
+
+ [A1,A2,Rho,Pi] = selectiontotale(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau);
+
+ %save output
+ dlmwrite(strcat(testFolder,'A1'), reshape(A1,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'A2'), reshape(A2,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'Rho'), reshape(Rho,1,[]), delimiter);
+ dlmwrite(strcat(testFolder,'Pi'), reshape(Pi,1,[]), delimiter);
+
+end
--- /dev/null
+function[] = testConstructionModelesLassoMLE()
+
+ testFolder = 'data/';
+ delimiter = '\n';
+
+ %get dimensions
+ dimensions = dlmread(strcat(testFolder,'dimensions'), delimiter);
+ n = dimensions(1);
+ p = dimensions(2);
+ m = dimensions(3);
+ k = dimensions(4);
+ L = dimensions(5);
+
+ %get all input arrays
+ phiInit = reshape(dlmread(strcat(testFolder,'phiInit'), delimiter), p, m, k);
+ rhoInit = reshape(dlmread(strcat(testFolder,'rhoInit'), delimiter), m, m, k);
+ piInit = transpose(dlmread(strcat(testFolder,'piInit'), delimiter));
+ gamInit = reshape(dlmread(strcat(testFolder,'gamInit'), delimiter), n, k);
+ mini = int64(dlmread(strcat(testFolder,'mini'), delimiter));
+ maxi = int64(dlmread(strcat(testFolder,'maxi'), delimiter));
+ gamma = dlmread(strcat(testFolder,'gamma'), delimiter);
+ glambda = dlmread(strcat(testFolder,'glambda'), delimiter);
+ X = reshape(dlmread(strcat(testFolder,'X'), delimiter), n, p);
+ Y = reshape(dlmread(strcat(testFolder,'Y'), delimiter), n, m);
+ seuil = dlmread(strcat(testFolder,'seuil'), delimiter);
+ tau = dlmread(strcat(testFolder,'tau'), delimiter);
+ A1 = int64(reshape(dlmread(strcat(testFolder,'A1'), delimiter), p, m+1, L));
+ A2 = int64(reshape(dlmread(strcat(testFolder,'A2'), delimiter), p, m+1, L));
+
+ %run constructionModelesLassoMLE.m
+ [phi,rho,pi,lvraisemblance] = constructionModelesLassoMLE(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2);
+
+ %get all stored outputs
+ ref_phi = reshape(dlmread(strcat(testFolder,'phi'), delimiter), p, m, k, L);
+ ref_rho = reshape(dlmread(strcat(testFolder,'rho'), delimiter), m, m, k, L);
+ ref_pi = reshape(dlmread(strcat(testFolder,'pi'), delimiter), k, L);
+ ref_lvraisemblance = reshape(dlmread(strcat(testFolder,'lvraisemblance'), delimiter), L, 2);
+
+ %check that output correspond to stored output
+ tol = 1e-5;
+ checkOutput('phi',phi,ref_phi,tol);
+ checkOutput('rho',rho,ref_rho,tol);
+ checkOutput('pi',pi,ref_pi,tol);
+ checkOutput('lvraisemblance',lvraisemblance,ref_lvraisemblance,tol);
+
+end
--- /dev/null
+function[] = testEMGLLF()
+
+ testFolder = 'data/';
+ delimiter = '\n';
+
+ %get dimensions
+ dimensions = dlmread(strcat(testFolder,'dimensions'), delimiter);
+ n = dimensions(1);
+ p = dimensions(2);
+ m = dimensions(3);
+ k = dimensions(4);
+
+ %get all input arrays
+ phiInit = reshape(dlmread(strcat(testFolder,'phiInit'), delimiter), p, m, k);
+ rhoInit = reshape(dlmread(strcat(testFolder,'rhoInit'), delimiter), m, m, k);
+ piInit = transpose(dlmread(strcat(testFolder,'piInit'), delimiter));
+ gamInit = reshape(dlmread(strcat(testFolder,'gamInit'), delimiter), n, k);
+ mini = int64(dlmread(strcat(testFolder,'mini'), delimiter));
+ maxi = int64(dlmread(strcat(testFolder,'maxi'), delimiter));
+ gamma = dlmread(strcat(testFolder,'gamma'), delimiter);
+ lambda = dlmread(strcat(testFolder,'lambda'), delimiter);
+ X = reshape(dlmread(strcat(testFolder,'X'), delimiter), n, p);
+ Y = reshape(dlmread(strcat(testFolder,'Y'), delimiter), n, m);
+ tau = dlmread(strcat(testFolder,'tau'), delimiter);
+
+ %run EMGLLF.m
+ [phi,rho,pi,LLF,S] = EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau);
+
+ %get all stored outputs
+ ref_phi = reshape(dlmread(strcat(testFolder,'phi'), delimiter), p, m, k);
+ ref_rho = reshape(dlmread(strcat(testFolder,'rho'), delimiter), m, m, k);
+ ref_pi = dlmread(strcat(testFolder,'pi'), delimiter);
+ ref_LLF = dlmread(strcat(testFolder,'LLF'), delimiter);
+ ref_S = reshape(dlmread(strcat(testFolder,'S'), delimiter), p, m, k);
+
+ %check that output correspond to stored output
+ tol = 1e-5;
+ checkOutput('phi',phi,ref_phi,tol);
+ checkOutput('rho',rho,ref_rho,tol);
+ checkOutput('pi',pi,ref_pi,tol);
+ checkOutput('LLF',LLF,ref_LLF,tol);
+ checkOutput('S',S,ref_S,tol);
+
+end
--- /dev/null
+function[] = testSelectiontotale()
+
+ testFolder = 'data/';
+ delimiter = '\n';
+
+ %get dimensions
+ dimensions = dlmread(strcat(testFolder,'dimensions'), delimiter);
+ n = dimensions(1);
+ p = dimensions(2);
+ m = dimensions(3);
+ k = dimensions(4);
+ L = dimensions(5);
+
+ %get all input arrays
+ phiInit = reshape(dlmread(strcat(testFolder,'phiInit'), delimiter), p, m, k);
+ rhoInit = reshape(dlmread(strcat(testFolder,'rhoInit'), delimiter), m, m, k);
+ piInit = transpose(dlmread(strcat(testFolder,'piInit'), delimiter));
+ gamInit = reshape(dlmread(strcat(testFolder,'gamInit'), delimiter), n, k);
+ mini = int64(dlmread(strcat(testFolder,'mini'), delimiter));
+ maxi = int64(dlmread(strcat(testFolder,'maxi'), delimiter));
+ gamma = dlmread(strcat(testFolder,'gamma'), delimiter);
+ glambda = dlmread(strcat(testFolder,'glambda'), delimiter);
+ X = reshape(dlmread(strcat(testFolder,'X'), delimiter), n, p);
+ Y = reshape(dlmread(strcat(testFolder,'Y'), delimiter), n, m);
+ seuil = dlmread(strcat(testFolder,'seuil'), delimiter);
+ tau = dlmread(strcat(testFolder,'tau'), delimiter);
+
+ %run constructionModelesLassoMLE.m
+ [A1,A2,Rho,Pi] = selectiontotale(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau);
+
+ %get all stored outputs
+ ref_A1 = int64(reshape(dlmread(strcat(testFolder,'A1'), delimiter), p, m+1, L));
+ ref_A2 = int64(reshape(dlmread(strcat(testFolder,'A2'), delimiter), p, m+1, L));
+ ref_Rho = reshape(dlmread(strcat(testFolder,'Rho'), delimiter), m, m, k, L);
+ ref_Pi = reshape(dlmread(strcat(testFolder,'Pi'), delimiter), k, L);
+
+ %check that output correspond to stored output
+ tol = 1e-5;
+ checkOutput('A1',A1,ref_A1,tol);
+ checkOutput('A2',A2,ref_A2,tol);
+ checkOutput('Rho',Rho,ref_Rho,tol);
+ checkOutput('Pi',Pi,ref_Pi,tol);
+
+end
--- /dev/null
+#include "constructionModelesLassoMLE.h"
+#include "ioutils.h"
+
+int main(int argc, char** argv)
+{
+ // read dimensions
+ const Int nbDims = 5;
+ Int* dimensions = readArray_int("dimensions",&nbDims,1);
+ mwSize n = dimensions[0];
+ mwSize p = dimensions[1];
+ mwSize m = dimensions[2];
+ mwSize k = dimensions[3];
+ mwSize L = dimensions[4];
+ free(dimensions);
+ mwSize lengthOne = 1;
+
+ ////////////
+ // INPUTS //
+ ////////////
+
+ // phiInit
+ const mwSize dimPhiInit[] = {p, m, k};
+ Real* phiInit = readArray_real("phiInit",dimPhiInit,3);
+
+ // rhoInit
+ const mwSize dimRhoInit[] = {m, m, k};
+ Real* rhoInit = readArray_real("rhoInit",dimRhoInit,3);
+
+ // piInit
+ Real* piInit = readArray_real("piInit",&k,1);
+
+ // gamInit
+ const mwSize dimGamInit[] = {n, k};
+ Real* gamInit = readArray_real("gamInit",dimGamInit,2);
+
+ // min number of iterations
+ Int* pmini = readArray_int("mini",&lengthOne,1);
+ Int mini = *pmini;
+ free(pmini);
+
+ // max number of iterations
+ Int* pmaxi = readArray_int("maxi",&lengthOne,1);
+ Int maxi = *pmaxi;
+ free(pmaxi);
+
+ // gamma
+ Real* pgamma = readArray_real("gamma",&lengthOne,1);
+ Real gamma = *pgamma;
+ free(pgamma);
+
+ // lambda
+ Real* glambda = readArray_real("glambda",&L,1);
+
+ // X
+ const mwSize dimX[] = {n, p};
+ Real* X = readArray_real("X",dimX,2);
+
+ // Y
+ const mwSize dimY[] = {n, m};
+ Real* Y = readArray_real("Y",dimY,2);
+
+ // seuil
+ Real* pseuil = readArray_real("seuil",&lengthOne,1);
+ Real seuil = *pseuil;
+ free(pseuil);
+
+ // tau
+ Real* ptau = readArray_real("tau",&lengthOne,1);
+ Real tau = *ptau;
+ free(ptau);
+
+ // A1
+ const mwSize dimA[] = {p, m+1, L};
+ Int* A1 = readArray_int("A1",dimA,3);
+
+ // A2
+ Int* A2 = readArray_int("A2",dimA,3);
+
+ /////////////
+ // OUTPUTS //
+ /////////////
+
+ // phi
+ const mwSize dimPhi[] = {dimPhiInit[0], dimPhiInit[1], dimPhiInit[2], L};
+ Real* phi = (Real*)malloc(dimPhi[0]*dimPhi[1]*dimPhi[2]*dimPhi[3]*sizeof(Real));
+
+ // rho
+ const mwSize dimRho[] = {dimRhoInit[0], dimRhoInit[1], dimRhoInit[2], L};
+ Real* rho = (Real*)malloc(dimRho[0]*dimRho[1]*dimRho[2]*dimRho[3]*sizeof(Real));
+
+ // pi
+ const mwSize dimPi[] = {k, L};
+ Real* pi = (Real*)malloc(dimPi[0]*dimPi[1]*sizeof(Real));
+
+ // lvraisemblance
+ const mwSize dimLvraisemblance[] = {L, 2};
+ Real* lvraisemblance = (Real*)malloc(dimLvraisemblance[0]*dimLvraisemblance[1]*sizeof(Real));
+
+ /////////////////////////////////////////
+ // Call to constructionModelesLassoMLE //
+ /////////////////////////////////////////
+
+ constructionModelesLassoMLE(
+ phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2,
+ phi,rho,pi,lvraisemblance,
+ n,p,m,k,L);
+
+ free(phiInit);
+ free(rhoInit);
+ free(piInit);
+ free(gamInit);
+ free(X);
+ free(Y);
+ free(A1);
+ free(A2);
+ free(glambda);
+
+ // Compare to reference outputs
+ Real* ref_phi = readArray_real("phi",dimPhi,4);
+ compareArray_real("phi", phi, ref_phi, dimPhi[0]*dimPhi[1]*dimPhi[2]*dimPhi[3]);
+ free(phi);
+ free(ref_phi);
+
+ // rho
+ Real* ref_rho = readArray_real("rho",dimRho,4);
+ compareArray_real("rho", rho, ref_rho, dimRho[0]*dimRho[1]*dimRho[2]*dimRho[3]);
+ free(rho);
+ free(ref_rho);
+
+ // pi
+ Real* ref_pi = readArray_real("pi",dimPi,2);
+ compareArray_real("pi", pi, ref_pi, dimPi[0]*dimPi[1]);
+ free(pi);
+ free(ref_pi);
+
+ // lvraisemblance
+ Real* ref_lvraisemblance = readArray_real("lvraisemblance",dimLvraisemblance,2);
+ compareArray_real("lvraisemblance", lvraisemblance, ref_lvraisemblance, dimLvraisemblance[0]*dimLvraisemblance[1]);
+ free(lvraisemblance);
+ free(ref_lvraisemblance);
+
+ return 0;
+}
--- /dev/null
+#include "EMGLLF.h"
+#include "ioutils.h"
+
+int main(int argc, char** argv)
+{
+ ////////////
+ // INPUTS //
+ ////////////
+
+ Int* dimensions = readArray_int("dimensions");
+ mwSize n = dimensions[0];
+ mwSize p = dimensions[1];
+ mwSize m = dimensions[2];
+ mwSize k = dimensions[3];
+ free(dimensions);
+
+ Real* phiInit = readArray_real("phiInit");
+ Real* rhoInit = readArray_real("rhoInit");
+ Real* piInit = readArray_real("piInit");
+ Real* gamInit = readArray_real("gamInit");
+ Int mini = read_int("mini");
+ Int maxi = read_int("maxi");
+ Real gamma = read_real("gamma");
+ Real lambda = read_real("lambda");
+ Real* X = readArray_real("X");
+ Real* Y = readArray_real("Y");
+ Real tau = read_real("tau");
+
+ /////////////
+ // OUTPUTS //
+ /////////////
+
+ Real* phi = (Real*)malloc(p*m*k*sizeof(Real));
+ Real* rho = (Real*)malloc(m*m*k*sizeof(Real));
+ Real* pi = (Real*)malloc(k*sizeof(Real));
+ Real* LLF = (Real*)malloc(maxi*sizeof(Real));
+ Real* S = (Real*)malloc(p*m*k*sizeof(Real));
+
+ ////////////////////
+ // Call to EMGLLF //
+ ////////////////////
+
+ EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau,
+ phi,rho,pi,LLF,S,
+ n,p,m,k);
+
+ free(phiInit);
+ free(rhoInit);
+ free(piInit);
+ free(gamInit);
+ free(X);
+ free(Y);
+
+ // Compare to reference outputs
+ Real* ref_phi = readArray_real("phi");
+ compareArray_real("phi", phi, ref_phi, p*m*k);
+ free(phi);
+ free(ref_phi);
+
+ Real* ref_rho = readArray_real("rho");
+ compareArray_real("rho", rho, ref_rho, m*m*k);
+ free(rho);
+ free(ref_rho);
+
+ Real* ref_pi = readArray_real("pi");
+ compareArray_real("pi", pi, ref_pi, k);
+ free(pi);
+ free(ref_pi);
+
+ Real* ref_LLF = readArray_real("LLF", maxi);
+ compareArray_real("LLF", LLF, ref_LLF);
+ free(LLF);
+ free(ref_LLF);
+
+ Real* ref_S = readArray_real("S");
+ compareArray_real("S", S, ref_S, p*m*k);
+ free(S);
+ free(ref_S);
+
+ return 0;
+}
--- /dev/null
+#include "EMGrank.h"
+#include "ioutils.h"
+
+int main(int argc, char** argv)
+{
+ // read dimensions
+ const Int nbDims = 4;
+ Int* dimensions = readArray_int("dimensions",&nbDims,1);
+ mwSize n = dimensions[0];
+ mwSize p = dimensions[1];
+ mwSize m = dimensions[2];
+ mwSize k = dimensions[3];
+ free(dimensions);
+ mwSize lengthOne = 1;
+
+ ////////////
+ // INPUTS //
+ ////////////
+
+ // Rho
+ const mwSize dimRho[] = {m, m, k};
+ Real* Rho = readArray_real("Rho",dimRho,3);
+
+ // Pi
+ Real* Pi = readArray_real("Pi",&k,1);
+
+ // min number of iterations
+ Int* pmini = readArray_int("mini",&lengthOne,1);
+ Int mini = *pmini;
+ free(pmini);
+
+ // max number of iterations
+ Int* pmaxi = readArray_int("maxi",&lengthOne,1);
+ Int maxi = *pmaxi;
+ free(pmaxi);
+
+ // X
+ const mwSize dimX[] = {n, p};
+ Real* X = readArray_real("X",dimX, 2);
+
+ // Y
+ const mwSize dimY[] = {n, m};
+ Real* Y = readArray_real("Y",dimY, 2);
+
+ // tau
+ Real* ptau = readArray_real("tau",&lengthOne,1);
+ Real tau = *ptau;
+ free(ptau);
+
+ // tau
+ Int* rank = readArray_int("rank",&k,1);
+
+ /////////////
+ // OUTPUTS //
+ /////////////
+
+ // phi
+ const mwSize dimPhi[] = {p, m, k};
+ Real* phi = (Real*)malloc(dimPhi[0]*dimPhi[1]*dimPhi[2]*sizeof(Real));
+
+ // LLF
+ Real* LLF = (Real*)malloc(1*sizeof(Real));
+
+ //////////////////////////
+ // Main call to EMGrank //
+ //////////////////////////
+
+ EMGrank(Pi,Rho,mini,maxi,X,Y,tau,rank,
+ phi,LLF,
+ n,p,m,k);
+
+ // free input pointers
+ free(Rho);
+ free(Pi);
+ free(X);
+ free(Y);
+ free(rank);
+
+ // Compare to reference outputs
+ Real* ref_phi = readArray_real("phi",dimPhi, 3);
+ compareArray_real("phi", phi, ref_phi, dimPhi[0]*dimPhi[1]*dimPhi[2]);
+ free(phi);
+ free(ref_phi);
+
+ // LLF
+ Real* ref_LLF = readArray_real("LLF",&lengthOne,1);
+ compareArray_real("LLF", LLF, ref_LLF, 1);
+ free(LLF);
+ free(ref_LLF);
+
+ return 0;
+}
--- /dev/null
+#include "constructionModelesLassoRank.h"
+#include "ioutils.h"
+
+int main(int argc, char** argv)
+{
+ // read dimensions
+ const Int nbDims = 5;
+ Int* dimensions = readArray_int("dimensions",&nbDims,1);
+ mwSize n = dimensions[0];
+ mwSize p = dimensions[1];
+ mwSize m = dimensions[2];
+ mwSize k = dimensions[3];
+ mwSize L = dimensions[4];
+ free(dimensions);
+ mwSize lengthOne = 1;
+
+ ////////////
+ // INPUTS //
+ ////////////
+
+ // piInit
+ const mwSize dimPi[] = {k, L};
+ Real* Pi = readArray_real("Pi",dimPi,2);
+
+ // rhoInit
+ const mwSize dimRho[] = {m, m, k, L};
+ Real* Rho = readArray_real("Rho",dimRho,4);
+
+ // min number of iterations
+ Int* pmini = readArray_int("mini",&lengthOne,1);
+ Int mini = *pmini;
+ free(pmini);
+
+ // max number of iterations
+ Int* pmaxi = readArray_int("maxi",&lengthOne,1);
+ Int maxi = *pmaxi;
+ free(pmaxi);
+
+ // X
+ const mwSize dimX[] = {n, p};
+ Real* X = readArray_real("X",dimX,2);
+
+ // Y
+ const mwSize dimY[] = {n, m};
+ Real* Y = readArray_real("Y",dimY,2);
+
+ // tau
+ Real* ptau = readArray_real("tau",&lengthOne,1);
+ Real tau = *ptau;
+ free(ptau);
+
+ // A1
+ const mwSize dimA[] = {p, L};
+ Int* A1 = readArray_int("A1",dimA,2);
+
+ // rangmin
+ Int* prangmin = readArray_int("rangmin",&lengthOne,1);
+ Int rangmin = *prangmin;
+ free(prangmin);
+
+ // rangmax
+ Int* prangmax = readArray_int("rangmax",&lengthOne,1);
+ Int rangmax = *prangmax;
+ free(prangmax);
+
+ /////////////
+ // OUTPUTS //
+ /////////////
+
+ // phi
+ mwSize Size = (mwSize)pow(rangmax-rangmin+1, k);
+ const mwSize dimPhi[] = {p, m, k, L*Size};
+ Real* phi = (Real*)malloc(dimPhi[0]*dimPhi[1]*dimPhi[2]*dimPhi[3]*sizeof(Real));
+
+ // lvraisemblance
+ const mwSize dimLvraisemblance[] = {L*Size, 2};
+ Real* lvraisemblance = (Real*)malloc(dimLvraisemblance[0]*dimLvraisemblance[1]*sizeof(Real));
+
+ //////////////////////////////////////////////
+ // Main call to constructionModelesLassoMLE //
+ //////////////////////////////////////////////
+
+ constructionModelesLassoRank(
+ Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax,
+ phi,lvraisemblance,
+ n,p,m,k,L);
+
+ free(Rho);
+ free(Pi);
+ free(X);
+ free(Y);
+ free(A1);
+
+ // Compare to reference outputs
+ Real* ref_phi = readArray_real("phi",dimPhi, 4);
+ compareArray_real("phi", phi, ref_phi, dimPhi[0]*dimPhi[1]*dimPhi[2]*dimPhi[3]);
+ free(phi);
+ free(ref_phi);
+
+ // lvraisemblance
+ Real* ref_lvraisemblance = readArray_real("lvraisemblance",dimLvraisemblance,2);
+ compareArray_real("lvraisemblance", lvraisemblance, ref_lvraisemblance, dimLvraisemblance[0]*dimLvraisemblance[1]);
+ free(lvraisemblance);
+ free(ref_lvraisemblance);
+
+ return 0;
+}
--- /dev/null
+#include "selectiontotale.h"
+#include "ioutils.h"
+
+mwSize main(mwSize argc, char** argv)
+{
+ // read dimensions
+ const Int nbDims = 5;
+ Int* dimensions = readArray_int("dimensions",&nbDims,1);
+ mwSize n = dimensions[0];
+ mwSize p = dimensions[1];
+ mwSize m = dimensions[2];
+ mwSize k = dimensions[3];
+ mwSize L = dimensions[4];
+ free(dimensions);
+ mwSize lengthOne = 1;
+
+ ////////////
+ // INPUTS //
+ ////////////
+
+ // phiInit
+ const mwSize dimPhiInit[] = {p, m, k};
+ Real* phiInit = readArray_real("phiInit",dimPhiInit,3);
+
+ // rhoInit
+ const mwSize dimRhoInit[] = {m, m, k};
+ Real* rhoInit = readArray_real("rhoInit",dimRhoInit,3);
+
+ // piInit
+ Real* piInit = readArray_real("piInit",&k,1);
+
+ // gamInit
+ const mwSize dimGamInit[] = {n, k};
+ Real* gamInit = readArray_real("gamInit",dimGamInit, 2);
+
+ // min number of iterations
+ Int* pmini = readArray_int("mini",&lengthOne,1);
+ Int mini = *pmini;
+ free(pmini);
+
+ // max number of iterations
+ Int* pmaxi = readArray_int("maxi",&lengthOne,1);
+ Int maxi = *pmaxi;
+ free(pmaxi);
+
+ // gamma
+ Real* pgamma = readArray_real("gamma",&lengthOne,1);
+ Real gamma = *pgamma;
+ free(pgamma);
+
+ // lambda
+ Real* glambda = readArray_real("glambda",&L,1);
+
+ // X
+ const mwSize dimX[] = {n, p};
+ Real* X = readArray_real("X",dimX,2);
+
+ // Y
+ const mwSize dimY[] = {n, m};
+ Real* Y = readArray_real("Y",dimY,2);
+
+ // seuil
+ Real* pseuil = readArray_real("seuil",&lengthOne,1);
+ Real seuil = *pseuil;
+ free(pseuil);
+
+ // tau
+ Real* ptau = readArray_real("tau",&lengthOne,1);
+ Real tau = *ptau;
+ free(ptau);
+
+ /////////////
+ // OUTPUTS //
+ /////////////
+
+ // A1
+ const mwSize dimA[] = {p, m+1, L};
+ Int* A1 = (Int*)malloc(dimA[0]*dimA[1]*dimA[2]*sizeof(Int));
+
+ // A2
+ Int* A2 = (Int*)malloc(dimA[0]*dimA[1]*dimA[2]*sizeof(Int));
+
+ // rho
+ const mwSize dimRho[] = {m, m, k, L};
+ Real* Rho = (Real*)malloc(dimRho[0]*dimRho[1]*dimRho[2]*dimRho[3]*sizeof(Real));
+
+ // pi
+ const mwSize dimPi[] = {k, L};
+ Real* Pi = (Real*)malloc(dimPi[0]*dimPi[1]*sizeof(Real));
+
+ //////////////////////////////////////////////
+ // Main call to constructionModelesLassoMLE //
+ //////////////////////////////////////////////
+
+ selectiontotale(
+ phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,
+ A1,A2,Rho,Pi,
+ n,p,m,k,L);
+
+ // free input pointers
+ free(phiInit);
+ free(rhoInit);
+ free(piInit);
+ free(gamInit);
+ free(glambda);
+ free(X);
+ free(Y);
+
+ // Compare to reference outputs
+ Int* ref_A1 = readArray_int("A1",dimA, 3);
+ compareArray_int("A1", A1, ref_A1, dimA[0]*dimA[1]*dimA[2]);
+ free(A1);
+ free(ref_A1);
+
+ // A2
+ Int* ref_A2 = readArray_int("A2",dimA, 3);
+ compareArray_int("A2", A2, ref_A2, dimA[0]*dimA[1]*dimA[2]);
+ free(A2);
+ free(ref_A2);
+
+ // Rho
+ Real* ref_Rho = readArray_real("Rho",dimRho, 4);
+ compareArray_real("Rho", Rho, ref_Rho, dimRho[0]*dimRho[1]*dimRho[2]*dimRho[3]);
+ free(Rho);
+ free(ref_Rho);
+
+ // Pi
+ Real* ref_Pi = readArray_real("Pi",dimPi, 2);
+ compareArray_real("Pi", Pi, ref_Pi, dimPi[0]*dimPi[1]);
+ free(Pi);
+ free(ref_Pi);
+
+ return 0;
+}
--- /dev/null
+// Check if array == refArray
+void compareArray(const char* ID, const void* array, const void* refArray, int size,
+ int isInteger)
+{
+ Real EPS = 1e-5; //precision
+ printf("Checking %s\n",ID);
+ Real maxError = 0.0;
+ for (mwSize i=0; i<size; i++)
+ {
+ Real error = isInteger
+ ? fabs(((Int*)array)[i] - ((Int*)refArray)[i])
+ : fabs(((Real*)array)[i] - ((Real*)refArray)[i]);
+ if (error >= maxError)
+ maxError = error;
+ }
+ if (maxError >= EPS)
+ printf(" Inaccuracy: max(abs(error)) = %g >= %g\n",maxError,EPS);
+ else
+ printf(" OK\n");
+}
+
+void compareArray_real(const char* ID, const void* array, const void* refArray, int size)
+{
+ return compareArray(ID, array, refArray, size, 0);
+}
+
+void compareArray_int(const char* ID, const void* array, const void* refArray, int size)
+{
+ return compareArray(ID, array, refArray, size, 1);
+}
+
+// Read array by columns (as in MATLAB) and return by-rows encoding
+void* readArray(const char* fileName, int isInteger)
+{
+ // need to prepend '../data/' (not really nice code...)
+ char* fullFileName = (char*)calloc(5+strlen(fileName)+1,sizeof(char));
+ strcat(fullFileName, "data/");
+ strcat(fullFileName, fileName);
+ FILE* file = fopen(fullFileName, "r");
+ free(fullFileName);
+
+ // first pass to know how many elements to allocate
+ // /////................... TODO
+
+ int d = 1;
+ for (int i=0; i<nbDims; i++)
+ totalDim *= dimensions[i];
+ size_t elementSize = isInteger
+ ? sizeof(Int)
+ : sizeof(Real);
+
+ // read all values, and convert them to by-rows matrices format
+ void* matlabArray = malloc(totalDim*elementSize);
+ char curChar = ' ';
+ char bufferNum[64];
+ for (mwSize u=0; u<totalDim; u++)
+ {
+ // position to next non-separator character
+ while (!feof(file) && (curChar==' ' || curChar=='\n' || curChar=='\t' || curChar==','))
+ curChar = fgetc(file);
+ // read number (as a string)
+ int bufferIndex = 0;
+ while (!feof(file) && curChar!=' ' && curChar!='\n' && curChar!='\t' && curChar!=',')
+ {
+ bufferNum[bufferIndex++] = curChar;
+ curChar = fgetc(file);
+ }
+ bufferNum[bufferIndex] = 0;
+ // transform string into Real, and store it at appropriate location
+ if (isInteger)
+ ((Int*)matlabArray)[u] = atoi(bufferNum);
+ else
+ ((Real*)matlabArray)[u] = atof(bufferNum);
+ }
+ fclose(file);
+
+ void* brArray = matlabToBrArray(matlabArray, dimensions, nbDims, isInteger);
+ free(matlabArray);
+ return brArray;
+}
+
projects / valse.git / commitdiff