gufodotto’s gists

gufodotto / server.R

Last active December 11, 2015 15:18

First Shiny webapp - Lorentz equation

	library(shiny)
	library(deSolve)
	library(ggplot2)

	# Define server logic required to plot various variables
	shinyServer(function(input, output) {

	solveLorenz <- function(pars, times=tout) {
	derivs <- function(t, state, pars) { # returns rate of change
	with(as.list(c(state, pars)), {

gufodotto / 6_Coll_Plot.r

Created July 8, 2012 08:02

	grid.newpage()
	# create a tile matrix with 0 and 1 depending on which parameter is being fitted.
	cl1<-ggplot(Coll_m) + geom_tile(aes(X1, X2, fill=value), , col='white') + scale_fill_gradient(low='white', high='black', limits=c(0, 1))
	# fill the same matrix with tiles of two different colors
	cl1<-cl1 + geom_text(aes(X1, X2, label=value, col=value)) + scale_color_gradient(low='black', high='white', limits=c(0, 1))
	cl1<-cl1 + opts(title="Parameters' N-ways Identifiability", axis.title.x = theme_blank(), axis.title.y = theme_blank(), axis.text.y = theme_text(colour = NA), axis.ticks = theme_segment(colour = NA), axis.text.x = theme_text(size = 12), legend.position='none')
	# now, on the right side, fill colors depending on the collinearity values
	cl2<-ggplot(Coll_f) + geom_tile(aes(X1, X2, fill=flag), col='white') + scale_fill_manual(values=c('green','red'))
	# finally print clipped collinearity values
	cl2<-cl2 + geom_text(data=Coll_f, aes(X1, X2, label=label), col='black') + theme_bw() + guides(fill=guide_legend(t

gufodotto / 5_MeltColl.r

Created July 8, 2012 08:02

	require(reshape)
	# select only the columns relative to parameter on/off, and the collinearity value
	Coll_m<-melt(t(Coll[,c(1:length(guess_pars),ncol(Coll))]));
	# split them in two different frames for easier handling
	Coll_f<-Coll_m[Coll_m$X1=='collinearity',];
	Coll_m<-Coll_m[Coll_m$X1!='collinearity',];
	# add a categorical label
	Coll_f$label[Coll_f$value<20]<- sprintf("%3.2g", Coll_f$value[Coll_f$value<20]); Coll_f$label[Coll_f$value>=20] <- ">= 20"; Coll_f$label[Coll_f$value>=100] <- "> 100";
	Coll_f$flag<-'ID OK'; Coll_f$flag[Coll_f$value>=20] <- "ID KO";
	Coll_f$flag <- factor(Coll_f$flag, levels=c('ID OK', 'ID KO'))

gufodotto / 4_Coll_out

Created July 8, 2012 08:01

gufodotto / 3_sensFun_n_Coll.r

Created July 8, 2012 08:00

	sf<-sensFun(func = Objective, parms = bullet_pars, varscale = 1)

	Coll <- collin(sF)

gufodotto / 2_plot_Lorentz.r

Created July 8, 2012 07:58

	pXY<-ggplot(as.data.frame(bullet)) +geom_path(aes(X, Y, alpha=Z), col='green') + opts(legend.position = "none")

	pXY<-pXY +geom_path(data=as.data.frame(target), aes(X, Y, alpha=Z), col='blue')

	pXY<-pXY +geom_path(data=as.data.frame(guess), aes(X, Y, alpha=Z), col='red')

	print(pXY)

gufodotto / 1_fit_guess_to_target.r

Created July 8, 2012 01:21

Lorentz_varia

	target_pars <- c(a = -9/3, b = -5, c = 30); target_pars
	target <- solveLorenz(target_pars)

	guess_pars<-c(a = -6/3, b = -8, c = 20); guess_pars
	guess <- solveLorenz(guess_pars)

	# print(system.time(Fit <- modFit(p = guess_pars, f = Objective))) # this works
	# print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="SANN", control=list(maxit=100), lower=c(-5,-10,0)))) # this works too
	# print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="SANN", control=list(maxit=10000)))) # this works too
	print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="Nelder-Mead", control=list(maxit=1000)))) # this works too

gufodotto / gist:3068849

Created July 8, 2012 01:10

Lorentz_ObjectiveFunction

	library(FME)

	Objective <- function(x, parset = names(x)) {
	guess_pars[parset] <- x
	tout <- seq(0, 50, by = 0.5)
	## output times
	bullet <- solveLorenz(guess_pars, tout)
	## Model cost
	return(modCost(obs = target, model = bullet))
	}

gufodotto / gist:3068827

Created July 8, 2012 00:57

Lorentz_short

	library(deSolve) # require this library

	solveLorenz <- function(pars, times=seq(0,100,by=0.1)) {
	derivs <- function(t, state, pars) { # returns rate of change
	with(as.list(c(state, pars)), {
	dX <- aX + YZ
	dY <- b * (Y-Z)
	dZ <- -XY + cY - Z
	return(list(c(dX, dY, dZ)))
	}

gufodotto / BubbleExoplanets.r

Created June 22, 2012 14:27

BubbleExoplanets.r

	# source("C:/Users/LucaF/Documents/My Dropbox/Software/R/Exoplanets.r")
	rm(list=ls(all=TRUE))

	if (!exists('palette_l')) {palette_l<-50} # how many steps in the palette

	N<-650 # number of exoplanets
	R<-1000 # radius of a 2D-circle where you'll want to plot them
	MinMass<-5; # minimum mass for the planets
	MaxMass<-100; # maximum mass for the planets
	P<-6 # power exponent for the mass distribution