Welcome, Guest

Author Topic: Coding  (Read 184419 times)

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #480 on: April 01, 2015, 12:42:31 PM »
tfw math homework is too much work. this only took me about twice the time it would take to do the homework by hand
---

an automated nonhomogenous linear ODE solver which uses the method of undetermined coefficients. the best part is that it outputs a latex file showing each step. includes batch processing

attached is the script, an example latex output, and the latex rendered as a pdf

Code: [Select]
from sympy import *
from copy import deepcopy

def all_terms(polyn, rcoeff = False):
    polyn = deepcopy(polyn)
    coeffs = []
    terms = []
    while polyn:
        fam, coeff = polyn.LT()
        famexpr = fam.as_expr()
        if famexpr == 1:
            fam = poly(x+1)-poly(x) #1
        else:
            fam = poly(fam.as_expr())
        terms.append(fam)
        coeffs.append(coeff)
        polyn = polyn - fam*coeff
    if rcoeff:
        return zip(terms, coeffs)
    else:
        return terms

def csub(pexp):
    pexp = deepcopy(pexp)
    for i in range(mc*2):
        pexp = pexp.subs('C'+str(i), 1)
    return pexp
       
def MOUC(problem, rlatex = True):
    latextext = []
    #print the problem
    latextext.append(latex(problem))
    latextext.append(latex("\\line(1,0){450}"))
   
    yc = Eq(problem.lhs, 0)
    latextext.append('y_c : ' + latex(yc))
    #print yc

    cp = yc
    for i, ndy in reversed(list(enumerate(dys))):
        cp = cp.subs(ndy, r**i)
    cp = cp.lhs
    latextext.append(latex(cp))

    cproots = []
    cpr = roots(Poly(cp))
    for k, v in cpr.items():
        cproots.append([k]*v)
    #print roots (imaginary!?)
    latextext.append('r : ' + ', '.join([latex(cproot) for cproot in cproots]))
    latextext.append(latex("\\line(1,0){450}"))
   
    ycsol = dsolve(yc.subs(y, y(x)))
    latextext.append(latex(ycsol))
    #print ycsol
   
    yfams = []
    ypf = poly(problem.rhs)
    yfams = all_terms(ypf)

    dyfams = [yfam.as_expr() for yfam in yfams]
    for yfam in dyfams:
        done = True
        while done:
            if type(yfam) in [type(S(val)) for val in xrange(-1,3)]:
                break
            dyfam = poly(diff(yfam, x))
            dyfammons = all_terms(dyfam)
            for dyfammon in dyfammons:
                dyfammon = dyfammon.as_expr()
                if dyfammon in dyfams:
                    done = False
                else:
                    dyfams.append(dyfammon)
                    dyfam = dyfammon
    dyfams1 = dyfams[:]
   
    yre = ycsol.rhs.expand()
    yre = csub(yre)
       
    rhsfams = [term.as_expr() for term in all_terms(poly(yre))]

    for i, dyfam in enumerate(dyfams):
        while dyfams[i] in rhsfams:
            dyfams[i] *= x
        while dyfams.count(dyfams[i]) > 1:
            dyfams[i] *= x
       
    #print promoted dyfams... if different
    for i, dyfam in enumerate(dyfams):
        if dyfam == dyfams1[i]:
            latextext.append(latex(dyfam))
        else:
            latextext.append(latex(dyfams1[i]) + ' \\rightarrow ' + latex(dyfam))

    Cs = [symbols('C'+str(i)) for i in range(3, 3+mc2)]
    yp = 0
    for i, term in enumerate(dyfams):
        yp += Cs[i]*term
       
    dyps = [yp]
    for i in xrange(1, md+1):
        dyps.append(diff(dyps[-1], x))
    #print each of the dyps
    latextext.append(latex("\\line(1,0){450}"))
    name = "y_p"
    for i in xrange(poly(cp).degree()+1):
        latextext.append(name + '=' + latex(dyps[i]))
        name += "'"

    adds = []
    TCs = all_terms(poly(cp), True)
    for T, C in TCs:
        adds.append(C*dyps[T.degree()])
    problemlhs = Add(*adds, evaluate=False)

    latextext.append(latex("\\line(1,0){450}"))
    ypbs = Eq(problemlhs, problem.rhs)
    ypas = Eq(problemlhs.collect(x), problem.rhs)
    #print unsimplified and simplified particular
    latextext.append(latex(ypbs))
    latextext.append(latex(ypas))
   
    eqs = {key: Eq(0,0, evaluate = False) for key in rhsfams+dyfams}
    for term, coeff in ypas.lhs.as_coefficients_dict().items():
        key = csub(term)
        cl, cr = eqs[key].lhs, eqs[key].rhs
        eqs[key] = Eq(cl+coeff*term, cr)
    for term, coeff in ypas.rhs.as_coefficients_dict().items():
        key = csub(term)
        cl, cr = eqs[key].lhs, eqs[key].rhs
        eqs[key] = Eq(cl, coeff*term+cr)
    for k, v in eqs.items():
        cl, cr = v.lhs, v.rhs
        eqs[k] = Eq((cl/k).simplify(), (cr/k).simplify())
    tosolve = []
    for k, v in eqs.items():
        if v is not S(True):
            tosolve.append(v)
           
    #print tosolve
    latextext.append(latex("\\line(1,0){450}"))
    for eq in tosolve:
        latextext.append(latex(eq))
           
    ucs = Cs[:len(tosolve)]
    Cps = {key:val for (key, val) in zip(ucs, solve_poly_system(tosolve, ucs)[0])}
    #print Cps
    for k, v in Cps.items():
        latextext.append(latex(k) + ' = ' + latex(v))
   
    ypsol = yp.xreplace(Cps)
    #print ypsol
    latextext.append(latex("\\line(1,0){450}"))
    latextext.append('y_p = ' +latex(ypsol))
   
    gsol = Eq(y, ycsol.rhs + ypsol)
    #print gsol
    latextext.append(latex(gsol))
   
    if rlatex:
        return gsol, latextext
    else:
        return gsol

def format(latexlines):
    start = "\\begin{dmath*}"
    end = "\\end{dmath*}"
    nl = "\n"
    return nl.join([start + nl + line + nl + end for line in latexlines])

def batchMOUC(problemdict, output):
    fulltext = ["\\documentclass[fleqn, 12pt]{article}\n\\usepackage{breqn}\n\\usepackage[cm]{fullpage}\n\\setlength{\\mathindent}{0pt}\n\\begin{document}"]
    #head, tail, and section
    for i, v in problemdict.items():
        _, latext = MOUC(v)
        fulltext.append('\\section*{'+str(i)+'}')
        fulltext.append(format(latext))
    fulltext.append("\\end{document}")
    with open(output, 'w') as ouf:
        ouf.write('\n'.join(fulltext))
       
init_printing()

md = mc = mc2 = 5
x, y, r = symbols('x y r')
dys = [y]
for i in xrange(1, md+1):
    dys.append(Derivative(dys[-1], x))
   
problems = {
    9: Eq(dys[2] + 2*dys[1] - 3*y, 1+x*exp(x)),
    10: Eq(dys[2] + 9*y, 2*cos(3*x) + 3*sin(3*x)),
    13: Eq(dys[2] + 2*dys[1] + 5*y, exp(x)*sin(x))
}

batchMOUC(problems, 'output.latex')

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #481 on: April 02, 2015, 07:48:21 AM »
number of "favorites" a piece gets on fanfiction.net vs which day it was published. i only considered stories that were published all at once (no incremental publishing). sample size for each day is about 90000. error bars are one standard deviation. monday is 0 and sunday is 6.

tldr publish on monday or thursday or friday. although it's less than a 2% difference from the average so it's probably not that important.

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #482 on: April 06, 2015, 11:01:01 AM »
python handles double negatives well


atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #483 on: April 14, 2015, 08:52:42 AM »
python 3.4 master race

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #484 on: April 14, 2015, 11:51:18 AM »
>3.4
boooo

FiahOwl

  • *****
  • Posts: 1234
  • This is, to give a dog and in recompense desire my dog again.
Re: Coding
« Reply #485 on: April 14, 2015, 11:59:26 AM »
#2.7

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #486 on: April 14, 2015, 02:54:41 PM »
3.4 handles unicode encoding automatically with no frustrations

which is what a programming language/compiler should do in 2015

it's the future

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #487 on: April 14, 2015, 03:24:57 PM »
but the 2.7 libraries!

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #488 on: April 19, 2015, 07:38:56 AM »
my IOint class is like my IOlist class. except this one doesn't use inheritance because i've been too lazy to get it working

Code: [Select]
class IOint():
    def __init__(self, filename):
        self.value = None
        self.filename = filename
    def set(self, val):
        self.value = val
        self.save()
    def save(self):
        with open(self.filename, 'wb') as ouf:
            ouf.write(str(self.value))
    def load(self):
        with open(self.filename, 'rb') as inf:
            self.value = int(inf.readline().strip())

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #489 on: April 27, 2015, 08:07:19 AM »
my western pacific statistical model...

it creates pretty believable tracks but the options are very limited (there's one track per square degree per month) and (probably) very innaccurate. the other issue is that many of the square degrees have no data and the storm will stop moving, so i need to create some interpolation method of some sort

it also only does tracks and doesn't forecast strength at all, but that is for a later addition

in the coming months i will hopefully be able to do some case studies

i could try to 'forecast' past storms too and see how well that goes
« Last Edit: April 27, 2015, 08:11:59 AM by atomic7732 »

Darvince

  • *****
  • Posts: 1837
  • 差不多
Re: Coding
« Reply #490 on: April 27, 2015, 08:09:58 AM »
i was gonna say you should be a meteorologist as a money making profession and then i realized

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #491 on: April 27, 2015, 08:42:32 AM »
http://imgur.com/a/uO6uV#0

this is the data that the program has so as you can see there are a lot of holes (except during peak season), but i'm not sure how i'd interpolate between them and determine where is interpolatable and where isn't

each character is how many storms have data for that square degree, and if it's greater than 9, A = 10, B = 11, R = 28, Y = 35 (yes there is one of those) etc
« Last Edit: April 27, 2015, 08:47:23 AM by atomic7732 »

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #492 on: April 27, 2015, 06:40:17 PM »
i'm studying for the computer science ap exam and google loves to troll me


atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #493 on: April 27, 2015, 07:55:09 PM »


« Last Edit: April 27, 2015, 09:20:28 PM by atomic7732 »

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #494 on: May 01, 2015, 08:42:23 AM »
http://stackoverflow.com/questions/24978052/interpolation-over-regular-grid-in-python

this appears to be very useful for filling in my grid but i have no idea how to use numpy/scipy grids and how to convert the data i already have (which is basically just a 3d list of list of lists) into a grid that can be interpolated by these modules

it doesn't help that the example has everything named with seemingly random letters!

Edit: potentially useful for actually you know getting numpy and scipy on python 3.4 http://continuum.io/downloads#py34
« Last Edit: May 01, 2015, 09:32:24 AM by atomic7732 »

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #495 on: May 01, 2015, 02:01:19 PM »
it's actually  not that hard. i never bothered with tutorials

converting:
np.array(yourlist) works

everything else is pretty much the same, except you can't append and can only store one type (like arrays are supposed to work).

slicing is all the same. iteration is all the same, BUT DON'T ITERATE
if you want to take the square of everything it's as simple as arr**2
or if you want to apply function f to everything it's as simple as arr = np.vectorize(f)(arr)

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #496 on: May 01, 2015, 02:04:41 PM »
kolok will try it later tonight

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #497 on: May 04, 2015, 08:02:19 AM »
ok so um this is not interpolation

This is what the original data looks like. It's visibly upside-down but that's fine, the program handles it properly.



And this is the "interpolated" result:



i somehow managed to flip the grid dimensions which i think reorganized the data and then didn't even interpolate anything



i fixed that and changed all the zero's to nan's because they were not being interpolated at all
and now i'm getting an error at
Code: [Select]
f = interpolate.Rbf(rr, cc, ar, function='linear')

ValueError: array must not contain infs or NaNs



yeeeeeeeeeeeeee



no time in class rn but this is the PTPC model's Tropical Storm Noul prediction

Code: [Select]
TS NOUL
[+0hr] +9.5 +139.6
[+6hr] +10.2 +140.5
[+12hr] +11.4 +140.7
[+18hr] +11.6 +140.9
[+24hr] +11.8 +141.1
[+30hr] +12.3 +141.4
[+36hr] +12.5 +141.7
[+42hr] +12.8 +141.9
[+48hr] +13.1 +142.5
[+54hr] +13.7 +143.1
[+60hr] +14.4 +143.9
[+66hr] +15.2 +144.7
[+72hr] +16.0 +145.6
[+78hr] +17.0 +146.6
[+84hr] +17.9 +147.8
[+90hr] +19.1 +149.2
[+96hr] +20.6 +150.7
[+102hr] +22.1 +152.5
[+108hr] +23.7 +154.4
[+114hr] +25.3 +156.0
[+120hr] +26.9 +156.7
[+126hr] +28.0 +157.0
[+132hr] +28.6 +158.2
[+138hr] +29.9 +159.5
[+144hr] +31.0 +160.0
[+150hr] +31.5 +161.2
[+156hr] +32.0 +161.7
[+162hr] +32.5 +163.2
[+168hr] +33.5 +164.5
[+174hr] +34.6 +165.7
[+180hr] +35.7 +166.7
forecast end
« Last Edit: May 04, 2015, 08:50:23 AM by atomic7732 »

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #498 on: May 04, 2015, 03:48:32 PM »
heeeeyyyyyy

the first functional prediction of an active storm

to do:
- inertial assumption for the first ~24 hrs, so as to preserve initial motion of the storm and hopefully put it on a more realistic track before it becomes pure statistics
- further division of the grid data into three types for El Niño, La Niña, and neutral years, which can affect storm track and intensity - though i'm not sure if there'd be too sparse of data, it might be a feature that only works during peak season (~Jul-Nov)
- intensity prediction (statistical)
- sst data which affects intensity, possible, maybe use anomalies, since sst averages are likely to be reflected in the overall storm data anyway
- shear data which affects intensity, unlikely to be implemented

blotz

  • Formerly 'bong'
  • *****
  • Posts: 813
  • op pls
Re: Coding
« Reply #499 on: May 04, 2015, 03:57:57 PM »
why is it on darrinward.com not kal's website

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #500 on: May 04, 2015, 04:08:38 PM »
because that's just an easy online lat-lon plotter that generates a link every time you use it

i could maybe later code a track plotter in matplotlib? maybe

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #501 on: May 05, 2015, 07:58:43 AM »
coding this makes me hot

hhnnnnggggg cyclone data

lol omg i added intensity prediction to my model and uh
it thinks Noul is going to be a Haiyan

Code: [Select]
TROPICAL STORM NOUL - 05 05 12Z
[+0hr] +9.5 +138.9 45.0 kt 992.0 hPa
[+6hr] +9.5 +138.9 53.0 kt 991.0 hPa
[+12hr] +9.7 +138.6 60.0 kt 989.0 hPa
[+18hr] +9.9 +138.0 68.0 kt 988.0 hPa
[+24hr] +10.3 +137.3 76.0 kt 987.0 hPa
[+30hr] +10.9 +136.3 79.0 kt 985.0 hPa
[+36hr] +11.4 +134.3 80.0 kt 983.0 hPa
[+42hr] +11.6 +133.2 84.0 kt 980.0 hPa
[+48hr] +11.9 +132.6 87.0 kt 980.0 hPa
[+54hr] +12.4 +132.2 89.0 kt 974.0 hPa
[+60hr] +12.8 +132.2 98.0 kt 971.0 hPa
[+66hr] +13.2 +132.1 106.0 kt 968.0 hPa
[+72hr] +13.9 +132.0 111.0 kt 962.0 hPa
[+78hr] +14.6 +132.0 116.0 kt 956.0 hPa
[+84hr] +15.0 +131.7 118.0 kt 952.0 hPa
[+90hr] +15.4 +131.5 119.0 kt 949.0 hPa
[+96hr] +16.4 +131.4 125.0 kt 941.0 hPa
[+102hr]+16.8 +131.2 130.0 kt 933.0 hPa
[+108hr]+17.1 +130.9 134.0 kt 926.0 hPa
[+114hr]+17.9 +130.3 139.0 kt 921.0 hPa
[+120hr]+18.7 +129.7 144.0 kt 916.0 hPa
[+126hr]+19.2 +129.3 154.0 kt 906.0 hPa
[+132hr]+19.9 +128.7 157.0 kt 902.0 hPa
[+138hr]+20.7 +128.5 162.0 kt 895.0 hPa
[+144hr]+21.6 +129.3 162.0 kt 897.0 hPa
[+150hr]+22.4 +130.8 159.0 kt 905.0 hPa
[+156hr]+23.4 +131.6 153.0 kt 910.0 hPa
[+162hr]+24.6 +132.7 150.0 kt 915.0 hPa
[+168hr]+25.7 +134.1 148.0 kt 918.0 hPa
[+174hr]+27.3 +135.6 143.0 kt 924.0 hPa
[+180hr]+28.5 +137.5 138.0 kt 927.0 hPa
forecast end

I feel like it's going to need to be modified, maybe to intensify only like 50 - 75% of what it is told to do but I'm not sure how to empirically decide this.

if i get the shear data working i could maybe have that "nerf" the intensity
« Last Edit: May 05, 2015, 08:25:29 AM by atomic7732 »

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #502 on: May 12, 2015, 08:31:46 AM »
yee beautiful

Code: [Select]
# http://www.darrinward.com/lat-long/?id=552821
# for plotting
 
# http://stackoverflow.com/questions/24978052/interpolation-over-regular-grid-in-python
# for scipy interpolation help
 
# https://knowledge.safe.com/articles/How_To/Calculating-accurate-length-in-meters-for-lat-long-coordinate-systems
# for lat/lon conversions to nautical miles
 
# http://matplotlib.org/basemap/users/examples.html
 
# create a 1 deg x 1 deg resolution grid for each month
# put SST data for each month
# put shear data?
 
#split SSTs by el nino and la nina and neutral
 
import math, sys
import numpy as np
import scipy.interpolate as interpolate
import matplotlib.pyplot as plt
 
class cycdat:
        def __init__(self, vx, vy, wind, pres, data):
                self.x = vx
                self.y = vy
                self.w = wind
                self.p = pres
                self.d = data #how many points of data are in this grid, so that velocities can be summed and averaged
                #self.t = temp
               
f = open('bst8089.txt', 'r+')
k = f.read()
 
lines = k.split('\n')
#ln1 = lines[1]
 
#make a grid or something
grid = []
for m in range(0, 12):
        temp2 = []
        for y in range(0, 50):
                temp = []
                for x in range(100, 180):
                        temp.append(cycdat(0, 0, 0, 0, 0))
                temp2.append(temp)
        grid.append(temp2)
       
#m-1, y, x-100
 
lastyear = None
 
for l in lines:
        sp = l.split(' ')
        sp = [y for y in sp if y != '']
       
        try:
                if sp[1] == '002':
                        #date
                        t1 = sp[0]
                        year1 = t1[:2]
                        month1 = int(t1[2:4])
                        day1 = t1[4:6]
         
                        #tropical or extratrop
                        if sp[2] == '6':
                                ex1 = False
                        else:
                                ex1 = True
               
                        #position
                        lat1 = int(sp[3])/10
                        lon1 = int(sp[4])/10
                 
                        #intensity
                        wind1 = int(sp[6])
                        pres1 = int(sp[5])
                       
                        if lastyear is not None and lastlon < 180.0 and lastlat < 50.0:
                                dlat = lat1 - lastlat
                                dlon = lon1 - lastlon
                                dwind = wind1 - lastwind
                                dpres = pres1 - lastpres
                               
                                grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100].x += dlat
                                grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100].y += dlon
                                grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100].w += dwind
                                grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100].p += dpres
                                grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100].d += 1
                       
                        lastyear = year1
                        lastmonth = month1
                        lastday = day1
                        lastlat = lat1
                        lastlon = lon1
                        lastwind = wind1
                        lastpres = pres1
                if sp[0] == '66666': #reset all of the last things and bypass adding a point between storm end & start positions
                        lastyear = None
        except:
                print("end of file... probably")
 
for list in grid:
        for list2 in reversed(list):
                for x in list2:
                        if x.d != 0:
                                n = x.d
                                if x.d > 9:
                                        n = chr(x.d + 55)
                                sys.stdout.write(str(n))
                                x.x = x.x/x.d
                                x.y = x.y/x.d
                                x.w = x.w/x.d
                                x.p = x.p/x.d
                        else:
                                sys.stdout.write(' ')
                sys.stdout.write('\n')
        sys.stdout.write('\n\n')
 
stormmonth = int(input('Initialization month: '))
stormlat = float(input('Initial latitude: '))
stormlon = float(input('Initial longitude: '))
stormwin = float(input('Initial wind speed (kt): '))
stormpre = float(input('Initial pressure (hPa): '))
stormsp = float(input('Initial speed (kt): ')) * 6 #6 hours of movement
stormdir = input('Initial direction (or STATIONARY): ')
#calculate nm/deg lon
rlat = stormlat * (math.pi/180)
hfactor = (111412.84 * math.cos(rlat) - 93.5 * math.cos(3*rlat))/1852
#calculate nm/deg lat
#calculate vector of speed/dir
if stormdir == "N":
    initv = stormsp / 60
    inith = 0
elif stormdir == "S":
    initv = -stormsp / 60
    inith = 0
elif stormdir == "W":
    inith = -stormsp / hfactor
    initv = 0
elif stormdir == "E":
    inith = stormsp / hfactor
    initv = 0
elif stormdir == "NE":
    inith = math.cos(math.pi/4) * stormsp / hfactor
    initv = math.sin(math.pi/4) * stormsp / 60
elif stormdir == "NW":
    inith = math.cos(3*math.pi/4) * stormsp / hfactor
    initv = math.sin(3*math.pi/4) * stormsp / 60
elif stormdir == "SW":
    inith = math.cos(5*math.pi/4) * stormsp / hfactor
    initv = math.sin(5*math.pi/4) * stormsp / 60
elif stormdir == "SE":
    inith = math.cos(7*math.pi/4) * stormsp / hfactor
    initv = math.sin(7*math.pi/4) * stormsp / 60
elif stormdir == "ENE":
    inith = math.cos(math.pi/8) * stormsp / hfactor
    initv = math.sin(math.pi/8) * stormsp / 60
elif stormdir == "NNE":
    inith = math.cos(3*math.pi/8) * stormsp / hfactor
    initv = math.sin(3*math.pi/8) * stormsp / 60
elif stormdir == "NNW":
    inith = math.cos(5*math.pi/8) * stormsp / hfactor
    initv = math.sin(5*math.pi/8) * stormsp / 60
elif stormdir == "WNW":
    inith = math.cos(7*math.pi/8) * stormsp / hfactor
    initv = math.sin(7*math.pi/8) * stormsp / 60
elif stormdir == "WSW":
    inith = math.cos(9*math.pi/8) * stormsp / hfactor
    initv = math.sin(9*math.pi/8) * stormsp / 60
elif stormdir == "SSW":
    inith = math.cos(11*math.pi/8) * stormsp / hfactor
    initv = math.sin(11*math.pi/8) * stormsp / 60
elif stormdir == "SSE":
    inith = math.cos(13*math.pi/8) * stormsp / hfactor
    initv = math.sin(13*math.pi/8) * stormsp / 60
elif stormdir == "ESE":
    inith = math.cos(15*math.pi/8) * stormsp / hfactor
    initv = math.sin(15*math.pi/8) * stormsp / 60
elif stormdir == "STATIONARY":
    inith = 0
    initv = 0
 
print('\n')
 
fi = open('output.txt', 'a')
 
hours = 0
print("[+" + str(hours).zfill(3) + "hr]\t+" + str(stormlat).zfill(4) + "\t+" + str(stormlon).zfill(5)  + "\t" + str(round(stormwin, 1)).zfill(5) + " kt\t" + str(round(stormpre, 1)).zfill(5) + " hPa")
fi.write(str(round(stormlat, 1)) + ", " + str(round(stormlon, 1)) + "\n")
 
hgrid = [] #horizontal (W-E) motion
vgrid = [] #vertical (N-S) motion
wgrid = [] #wind
pgrid = [] #pressure
for list in grid:
    htemp = []
    vtemp = []
    wtemp = []
    ptemp = []
    for list2 in list:
        htemp2 = []
        vtemp2 = []
        wtemp2 = []
        ptemp2 = []
        for x in list2:
            htemp2.append(x.x)
            vtemp2.append(x.y)
            wtemp2.append(x.w)
            ptemp2.append(x.p)
        htemp.append(htemp2)
        vtemp.append(vtemp2)
        wtemp.append(wtemp2)
        ptemp.append(ptemp2)
    hgrid.append(htemp)
    vgrid.append(vtemp)
    wgrid.append(wtemp)
    pgrid.append(ptemp)
   
def interp(g):
    ar = np.array(g[stormmonth - 1])
    ar[ ar==0 ] = np.nan
    r = np.linspace(0, 1, ar.shape[1])
    c = np.linspace(0, 1, ar.shape[0])
     
    rr, cc = np.meshgrid(r, c)
    vals = ~np.isnan(ar)
    f = interpolate.Rbf(rr[vals], cc[vals], ar[vals], function='linear')
    return(f(rr, cc))
   
hi = interp(hgrid) #INTERPOLATE HORIZONTAL GRID
vi = interp(vgrid) #INTERPOLATE VERTICAL GRID
wi = interp(wgrid) #INTERPOLATE WIND SPEED GRID
pi = interp(pgrid) #INTERPOLATE PRESSURE GRID
 
#plt.imshow(wi)
#plt.show()
 
def simulate(slat, slon, g, g2, g3, g4):
        tlat = g[math.floor(slat)][math.floor(slon)-100]
        tlon = g2[math.floor(slat)][math.floor(slon)-100]
        twin = g3[math.floor(slat)][math.floor(slon)-100]
        tpre = g4[math.floor(slat)][math.floor(slon)-100]
        return (tlat, tlon, twin, tpre)
 
for x in range(0, 180, 6):
        templat, templon, tempwind, temppres = simulate(stormlat, stormlon, hi, vi, wi, pi)
        if hours < 36:
            frac = (1/1296) * (hours - 36)**2
            stormlat += frac * initv + (1-frac) * templat
            stormlon += frac * inith + (1-frac) * templon
        else:
            stormlat += templat
            stormlon += templon
        stormwin += .5 * tempwind
        stormpre += .5 * temppres
        if stormlat == lastlat and stormlon == lastlon:
                print("storm stalled")
                break
        if stormlat > 50.0 or stormlon > 180.0 or stormlat < 0.0 or stormlon < 100.0:
                print("storm left grid")
                break
        lastlat = stormlat
        lastlon = stormlon
        hours += 6
        print("[+" + str(hours).zfill(3) + "hr]\t+" + str(round(stormlat, 1)).zfill(4) + "\t+" + str(round(stormlon, 1)).zfill(5) + "\t" + str(round(stormwin, 1)).zfill(5) + " kt\t" + str(round(stormpre, 1)).zfill(5) + " hPa")
        fi.write(str(round(stormlat, 1)) + ", " + str(round(stormlon, 1)) + "\n")
        if hours == 180:
                print("forecast end")
               
fi.close()

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #503 on: May 12, 2015, 11:33:40 AM »
i feel like there should be a better way to handle that block of elif statements for every direction you had

atomic7732

  • Global Moderator
  • *****
  • Posts: 3829
  • caught in the river turning blue
    • Paladin of Storms
Re: Coding
« Reply #504 on: May 12, 2015, 02:16:42 PM »
i agree but can't think of anything

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #505 on: May 12, 2015, 03:39:40 PM »
well i don't know how it works so i'm not going to try deciphering it but i think this piece of code should help.

Code: [Select]
def gettheta(tlc): #threelettercode
dd = {'E':0, 'N':1, 'W':2, 'S':3}
if 'S' in tlc and 'E' in tlc:
dd['E'] = 4
theta = dd[tlc[-1]]
for l in reversed(tlc):
theta = (dd[l] + theta)/2.
return theta/2.*math.pi #convert to radians

then you can do

Code: [Select]
inith = math.cos(gettheta(tlc)) * stormsp / hfactor
initv = math.sin(gettheta(tlc)) * stormsp / 60

or something

you still need an separate case for the STATIONARY

---

or for some shorter (but less elegant) code

Code: [Select]
dirmap = dict(zip(['E', 'ENE', 'NE', 'NNE', 'N', 'NNW', 'NW', 'WNW', 'W', 'WSW', 'SW', 'SSW', 'S', 'SSE', 'SE', 'ESE'], [i/8.*math.pi for i in range(16)]))
and

Code: [Select]
inith = math.cos(dirmap[tlc]) * stormsp / hfactor
initv = math.sin(dirmap[tlc]) * stormsp / 60

----

also considering that objects are passed by reference i'm pretty sure instead of

Code: [Select]
grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100].x += dlat
grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100].y += dlon
grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100].w += dwind
grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100].p += dpres
grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100].d += 1

you can do

Code: [Select]
cell = grid[lastmonth - 1][math.floor(lastlat)][math.floor(lastlon) - 100]
cell.x += dlat
cell.y += dlon
cell.w += dwind
cell.p += dpres
cell.d += 1
« Last Edit: May 12, 2015, 03:47:43 PM by vh »

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #506 on: May 25, 2015, 11:35:36 AM »
A tool I made for designing geometric minecraft builds. Plots ellipses, splines, parabolas, hyperbolas, and anything you can write as a parametric function.

Code: [Select]
def demo():
B = UnitGrid(-100, 100, -100, 100, False, False)
print 'plotting hyperbola'
B.hyperbola(17, 23, 54, -71, 1)
print 'plotting ellipse'
B.ellipse(-20, -80, 15, 17)
print 'plotting parabola'
B.parabola(0.01, -2, 20)
print 'plotting line'
B.line(54, -90, -90, 11)
print 'plotting spline'
B.spline([0,10,72,-4,-90,-100],[0,10,-11,17,-90,-100])
B.write('demo')

demo()

you can clearly see the ellipse, hyperbola, parabola, spline, and line in this image.

vh

  • formerly mudkipz
  • *****
  • Posts: 1138
  • "giving heat meaning"
Re: Coding
« Reply #507 on: June 01, 2015, 03:47:39 PM »
this is more stats than coding but whatever.

x axis: days since epoch
y axis: seconds since 0000
color: keypresses per five minutes

Darvince

  • *****
  • Posts: 1837
  • 差不多
Re: Coding
« Reply #508 on: June 01, 2015, 04:22:10 PM »
whenpulse

FiahOwl

  • *****
  • Posts: 1234
  • This is, to give a dog and in recompense desire my dog again.
Re: Coding
« Reply #509 on: June 01, 2015, 04:42:12 PM »
pyclick