Help:=proc()
print(`Lasso(data,lambda,error-tolerance) ExtractNice(arrInputVar,outputVar)`):


end:

#Data will be of form [[[a,b,c,1],y],[[a',b',c',1],y'],...]
Lasso:=proc(data,lam,er) local i,n,m,z,temp,w,wn,e,k,s,p,j:
n:=nops(data):
m:=nops(data[1][1]):
z:=[0$m]:
for j from 1 to m do
temp:=0:
for i from 1 to n do
temp:=temp+(data[i][1][j])^2:
od:
z[j]:=evalf(temp):
#z[j]:=1:
od:
w:=[0$m]:
wn:=[0$m]:
e:=er+1:

for k from 1 to 1000 while e>er do

s:=[0$n]:
for i from 1 to n do
temp:=0:
for j from 1 to m do 
temp:=temp+w[j]*data[i][1][j]:
od:
s[i]:=temp:
od:
for j from 1 to m do
temp:=0:
for i from 1 to n do
temp:=temp+data[i][1][j]*(data[i][2]-s[i]+w[j]*data[i][1][j]):
od:
#wn[j]:=temp/z[j]:
if temp<-1*lam/2 then 
wn[j]:=(temp+lam/2)/z[j]:
fi:
if temp >= -1*lam/2 and temp <= lam/2 then
wn[j]:=0:
fi:
if temp>lam/2 then
wn[j]:=(temp-lam/2)/z[j]:
fi:
od:
#print(evalf(wn)):
e:=add(seq((wn[j]-w[j])^2,j=1..m)):
w:=wn:
if k=10000 then 
print('bad'):
fi:
od:
wn:
end:


#####Here and below was just for testing things
RandData:=proc(n,m,e) local a,b,i,j,d,temp,ra,ra1:
a:=[0$(m)]:
ra:=rand(-10..10.0):
for i from 1 to m do 
a[i]:=ra():
od:
ra1:=rand(-e..e):
b:=[0$n]:
for i from 1 to n do
temp:=[0$(m-1),1]:
for j from 1 to m-1 do
temp[j]:=ra()/j:
od:
b[i]:=[temp,add(seq(temp[j]*a[j],j=1..m))+ra1()]:
od:
print(a):
b:
end:

MR:=
[[Balaban, Tadeusz, 3, 57, 240, 8, 136839,0, 1967], 
[Beals, Robert, 4, 30, 294, 10, 242262, 0, 1980], 
[Beck, Jozsef, 4, 120, 1043, 18, 190938,1, 1974], 
[Borisov, Lev, 3, 48, 688, 15, 138117,1,1996], 
[Buch, Anders, 3, 40, 632, 14, 127510,0,1999], 
[Cakoni, Fioralba, 3, 94, 1746, 24, 156319,1,1996], 
[Carbone, Lisa, 3, 30, 118, 6, 118040,0,1997], 
[Carlen, Eric, 4, 130, 1929, 23, 181415,1,1984], 
[Chanillo, Sagun, 4, 71, 1111, 18, 156419,1,1980], 
[Han, ZhengChao, 3, 19, 507, 11, 116846,0,1989], 
[Feehan, Paul, 3, 27, 198, 10, 142972,1,1992], 
[Gibney, Angela, 3, 16, 179,8, 142000,0,2000], 
[Goldin, Gerald, 4, 79, 193, 9, 200175,0,1969], 
[Goldstein, Sheldon, 4, 108, 641, 13, 162850,0,1974], 
[Huang, Xiaojun, 4, 66, 805, 17, 156388,1,1994], 
[Huang, YiZhi, 3, 62, 1617, 20, 135147,0,1990], 
[Iwaniec, Henryk, 4, 157, 4870, 32, 274968,1,1972], 
[Kahn, Jeffry, 4, 110, 1016, 16, 202922,1,1979], 
[Kiessling, Michael, 3, 62, 548, 11, 126634,1,1989], 
[Komlos, Janos, 4, 101, 2421, 27, 189158,1,1967], 
[Kontorovich, Alex, 3, 33, 228, 10, 170000,1,2007], 
[Kopparty, Swastik, 2, 58, 218, 8, 131432,0,2011], 
[Krashen, Daniel, 3, 21, 142, 7, 142000,1,2001],
[Kriventsov, Dennis, 1, 7, 38, 4, 100000,0,2015], 
[Larsen, Kasper, 2, 13, 86, 5, 177640,0,2004], 
[Lebowitz, Joel, 4, 380, 3493, 27, 274330,1,1956], 
[Lepowsky, James, 4, 96, 3562, 27, 198213,1, 1970], 
[Li, Yanyan, 4, 155, 4466, 35, 196762,1, 1988], 
[Liu, Liping, 2, 10, 24, 3, 132901,0, 2006], 
[Luo, Feng, 4, 88, 903, 14, 122705,1, 1989], 
[Miller, Stephen, 3, 37, 378, 12, 143601,1, 1997], 
[Mirek, Mariusz, 1, 18, 104, 6, 102000, 0, 2011], 
[Mischaikow, Konstantin, 4, 116, 2057, 25, 195136,1,1985], 
[Narayanan, Bhargav, 1, 20, 34, 3, 102000,0, 2015], 
[Nussbaum, Roger, 4, 138, 2530, 26, 185683,1,1969], 
[Ocone, Daniel, 3, 41, 462, 12, 139443,0,1980], 
[Retakh, Vladimir, 3, 85, 936, 13, 139525,1,1973], 
[Roberts, Fred, 4, 127, 963, 17, 242906,1, 1968], 
[Rong, Xiaochun, 4, 43, 340, 10, 146108,1,1990], 
[Sahi, Siddhartha,3, 61, 746, 15, 128875,0, 1985], 
[Saks, Michael, 4, 140, 1054, 14, 194051,0,1980], 
[Saraf, Shubhangi, 2, 45, 200, 9, 93240,0, 2012],
[Sargsyan, Grigor, 2, 26, 77, 5, 106329,0,2009], 
[Scheffer, Vladimir, 3, 15, 545, 9, 108125,0,1974], 
[Sesum, Natasa, 3, 40, 601, 15, 126919,1,2004], 
[Soffer, Avraham, 4, 81, 1461, 15, 156858,1,1984], 
[Song, Jian, 3, 37, 730, 16, 121329,0,2004], 
[Sun,Hongbin, 1, 15, 52, 5, 96000,0,2014], 
[Sussmann, Hector, 4, 119, 1755, 22, 186954,1,1969], 
[TahvildarZadeh, Shadi, 3, 25, 621, 12, 132894,0,1991], 
[Thomas, Simon, 4, 86, 557, 13, 207659,1,1983], 
[Tiep, PhamHuu, 4, 144, 1431, 19, 195000,1,1988], 
[Tunnell, Jerrold, 2, 9, 269, 6, 112012,1,1977], 
[Vogelius, Michael, 4, 80, 2561, 28, 288356,1,1980], 
[Weibel, Charles, 4, 115, 2893, 20, 179121,1,1977], 
[Weston, Kimberly, 1, 3, 3, 1, 100000,0,2016], 
[Woodward, Christopher, 4, 43, 807, 16, 153567,1,1996], 
[Zeilberger, Doron, 4, 215, 2432, 23, 217337,1,1976], 
[Zhang,Wujun, 1, 16, 69, 5, 94500,0,2012]]
:

#SubLists(L,I1,J1): given a list L of lists where in each item, the I1-th list describes a categorical featues from 1 to J1
#outputs the individual lists for each of the categorical features. Try:
#SubLists(MR,3,4);
SubLists:=proc(L,I1,J1) local T,J1A,K:

if not (type(L,list)  and {seq(type(L[K],list),K=1..nops(L))}={true} and nops({seq(nops(L[K]),K=1..nops(L))})=1 ) then
 print(`Bad input`):
 RETURN(FAIL):
fi:


if {seq(L[K][I1],K=1..nops(L))} minus {seq(K,K=1..J1)}<>{} then
 print(`Bad input`):
 RETURN(FAIL):
fi:
  

for J1A from 1 to J1 do
  T[J1A]:=[]:
od:

for K from 1 to nops(L) do
 T[L[K][I1]]:=  [op(T[L[K][I1]]),L[K]]:
od:

[seq(T[K],K=1..J1)]:
end:


#ExtractFields(L,M): inputs a list of lists L and a list M , produces a a list with only the features given by M in that order
#and sorted
#For example, try:
#ExtractFields(MR,[7,1,2]);
ExtractFields:=proc(L,M) local K1,K2:

sort([seq([seq(L[K1][M[K2]],K2=1..nops(M))],K1=1..nops(L))]):

end:

#AveAndSD(L): the average and s.d. of L
AveAndSD:=proc(L) local i,mu,sig:
mu:=evalf(add(L[i],i=1..nops(L))/nops(L)):
sig:=sqrt(add((L[i]-mu)^2,i=1..nops(L))/nops(L)):
[mu,sig]:
end:

#Cor(L): inputs a list of pairs [P,S] outputs the corrletion
Cor:=proc(L) local L1,L2,i,A,B,mu1,mu2,sig1,sig2:
L1:=[seq(L[i][1],i=1..nops(L))]:
L2:=[seq(L[i][2],i=1..nops(L))]:
A:=AveAndSD(L1):
B:=AveAndSD(L2):
mu1:=A[1]: sig1:=A[2]: mu2:=B[1]: sig2:=B[2]:

[A,B,add((L1[i]-mu1)*(L2[i]-mu2),i=1..nops(L))/(nops(L)*sig1*sig2)]:
end:

ExtractNice:=proc(L,y):
[seq([[seq(MR[K1][L[K2]],K2=1..nops(L))],MR[K1][y]],K1=1..nops(MR))]:
end:


ExtractNice2:=proc(L,y):
[seq([[seq(MR[K1][L[K2]],K2=1..nops(L)),1],MR[K1][y]],K1=1..nops(MR))]:
end:

RandData2:=proc(n,m,e) local a,b,i,j,d,temp,ra,ra1:
a:=[0$(m)]:
ra:=rand(-10..10.0):
for i from 1 to m do 
a[i]:=ra():
od:
ra1:=rand(-e..e):
b:=[0$n]:
for i from 1 to n do
temp:=[0$(m-1),1]:
for j from 1 to m-1 do
temp[j]:=ra()*j^3:
od:
b[i]:=[temp,add(seq(temp[j]*a[j],j=1..m))+ra1()]:
od:
print(a):
b:
end: