netfit_function(data,meshsizes,type="norm.loc",rel) {
 lens_rep(data[,1],ncol(data[,-1]))
 msizes_rep(meshsizes,rep(nrow(data),ncol(data[,-1])))
 dat_as.vector(data[,-1])
 var1_lens*msizes; var2_msizes^2; var3_(lens/msizes)^2 
 var4_lens/msizes; var5_log(msizes); var6_log(msizes/msizes[1])
 var7_var6*log(lens) - 0.5*var6*var6; var8_lens*lens ;
 var9_msizes/lens ;
 if(missing(rel))
  os_0
 else
  os_rep(log(rel),rep(nrow(data),ncol(data[,-1])))
 if(type=="norm.loc") {
  if(missing(rel))
   fit_glm(dat ~ var1 + var2 + as.factor(lens),family=poisson)
  else
   fit_glm(dat ~ var1 + var2 + as.factor(lens) + offset(os),family=poisson)
  x_fit$coef[c("var1","var2")] 
  pars_c(-msizes[1]*2*x[2]/x[1],sqrt(-2*x[2]/(x[1]*x[1])))   }
 if(type=="norm.sca") {
  if(missing(rel))
   fit_glm(dat ~ var3 + var4 + as.factor(lens),family=poisson)
  else
   fit_glm(dat ~ var3 + var4 + as.factor(lens) + offset(os),family=poisson)
  x_fit$coef[c("var3","var4")] 
  pars_c(-msizes[1]*x[2]/(2*x[1]),sqrt(-msizes[1]^2/(2*x[1])))       }
 if(type=="gamma") {
  if(missing(rel))
   fit_glm(dat ~ var4 + var5 + as.factor(lens),family=poisson)
  else
   fit_glm(dat ~ var4 + var5 + as.factor(lens) + offset(os),family=poisson)
  x_fit$coef[c("var4","var5")] 
  pars_c(-x[2]+1,-msizes[1]/x[1])}
 if(type=="lognorm") {
  if(missing(rel))
   fit_glm(dat ~ var6 + var7 + as.factor(lens),family=poisson)
  else
   fit_glm(dat ~ var6 + var7 + as.factor(lens) + offset(os),family=poisson)
  x_fit$coef[c("var6","var7")] 
  pars_c(-(x[1]-1)/x[2],sqrt(1/x[2])) }
 if(type=="inv.Gaussian"){
  if (missing(rel)){
   fit_glm(dat ~ 
      var4 + var9 +offset(3/2*log(msizes)) + as.factor(lens) ,family=poisson)}
  else{
   fit_glm(dat ~ 
      var4 + var9 +offset(3/2*log(msizes)+os) +as.factor(lens) ,family=poisson)}
  x_fit$coef[c("var4","var9")]
  pars_c(msizes[1]*sqrt(x[2]/x[1]),-2*msizes[1]*x[2])   }
 devres_round(matrix(resid(fit,type="deviance"),nrow(data),ncol(data[,-1])),4)
 g.o.f_c(fit$deviance,round(fit$df),fit$null)
 fit.type_paste(paste(type,ifelse(missing(rel),"",": with offset")))
 return(fitted.pars=x,fit.type=fit.type,x=pars,g.o.f,devres)        }


#Look at residuals matrix
plot.resids_
 function(residuals,msizes,lens,cex=1,title="Deviance residuals",...) {
 if(missing(lens)) lens_1:nrow(residuals)
 if(missing(msizes)) msizes_1:ncol(residuals)
 plot(1,1,xlim=c(min(lens),max(lens)),xlab="Length",ylab="Mesh size",
     ylim=c(min(msizes),max(msizes))+(cex/50)*c(-1,1)*(max(msizes)-min(msizes)),
     yaxp=c(min(msizes),max(msizes),length(msizes)-1),
     xlab="Length",ylab="Mesh size",type="n",main=title,...)
 for(i in 1:nrow(residuals))
  for(j in 1:ncol(residuals))
   points(lens[i],msizes[j],pch=ifelse(residuals[i,j]>0,16,1),
                            mkh=abs(residuals[i,j])*cex/100) }
  
#Example of use
#holt.dat_matrix(scan("holt.dat"),byrow=T,ncol=9)
#meshsizes_c(13.5,14.0,14.8,15.4,15.9,16.6,17.8,19.0)
#fit1_netfit(holt.dat,meshsizes,type="norm.loc")
#fit2_netfit(holt.dat,meshsizes,type="norm.loc",rel=meshsizes)
#X11(); par(mfrow=c(2,1))
#plot.resids(fit1$devres,meshsizes,holt.dat[,1],cex=3,title="")
#plot.resids(fit2$devres,meshsizes,holt.dat[,1],cex=3,title="")