gufodotto

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

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

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

> Coll
  a b c N collinearity
1 1 1 0 2          1.1
2 1 0 1 2          2.8
3 0 1 1 2          1.0
4 1 1 1 3          2.9
> 

gufodotto

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

Coll <- collin(sF)

gufodotto

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

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

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

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 <- a*X + Y*Z
				dY <- b * (Y-Z)
				dZ <- -X*Y + c*Y - Z
				return(list(c(dX, dY, dZ)))
			}

gufodotto

# 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